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Abstract 

We survey approaches to nonrelativistic density functional theory (DFT) for nuclei using 
progress toward ab initio DFT for Coulomb systems as a guide. Ab initio DFT starts with a micro- 
scopic Hamiltonian and is naturally formulated using orbital-based functionals, which generalize 
the conventional local-density-plus-gradients form. The orbitals satisfy single-particle equations 
with multiplicative (local) potentials. The DFT functionals can be developed starting from inter- 
nucleon forces using wave-function based methods or by Legendre transform via effective actions. 
We describe known and unresolved issues for applying these formulations to the nuclear many- 
body problem and discuss how ab initio approaches can help improve empirical energy density 
functionals. 

Keywords: Density functional theory, nuclear structure, many-body perturbation theory 



Contents 

1 Introduction 

1.1 Overview 

1.2 Basic features/ingredients of DFT 

1.3 Coulomb vs. nuclear DFT _ 

1.4 Scope and plan of review [l^ 

2 Orbital-based DFT 

2.1 Hartree-Fock in coordinate representation 

2.2 Motivation for orbital-dependent functionals 

2.3 Derivation of the optimized effective potential 

2.4 OEP from total energy minimization or density invariance 

2.5 Approximations 

3 DFT and ab initio wave function methods 

3.1 Goldstone many-body perturbation theory 

3.2 Improved perturbation theory 

3.3 Low- momentum interactions 

3.4 Density matrix expansion 

4 DFT as Legendre transform 

4.1 Analogy to Legendre transform in thermodynamics 

4.2 Effective actions for composite operators 

4.3 EFT and power counting for functionals 

4.4 Additional comments 



1 



5 Topics for nuclear DFT [49 



5.1 Pairing |49 



5.2 Broken symmetries |53 



5.3 Single-particle energies l55 



5.4 Improving empirical EDF's |56 



6 Summary and outlook |58 



References 62 



1 Introduction 
1.1 Overview 

Density functional theory (DFT) has been applied to the Coulomb many-body problem with great 
phenomenological success in predicting properties of atoms, molecules, and solids [H [21 [3l SI El [6] . DFT 
calculations are comparatively simple to implement yet often very accurate and have a computational 
cost that makes them at present the only choice for systems with large numbers of electrons [7] . For these 
same reasons (with nucleons rather than electrons), large-scale collaborations of nuclear physicists in the 
SciDAC UNEDF E] ("Universal Nuclear Energy Density Functional") and FIDIPRO pTOJ projects, 
as well as many other individuals, are working on further developing DFT for the nuclear many-body 
problem. Questions in astrophysics and the advent of new experimental facilities to study nuclei at 
the limits of existence, as well as societal needs, are driving multi-pronged efforts to calculate nuclear 
properties and reactions across the full table of the nuclides more accurately and reliably than what is 
currently possible with existing energy density functional (EDF) methods (e.g., those based on Skyrme, 
Gogny, or relativistic mean- field functionals [TT]). 

A principal strategy is to exploit the substantial and ongoing progress in ab initio nuclear struc- 
ture calculations, which are primarily based on approximating the many-nucleon wave function. This 
progress is the consequence of synergistic advances in the construction of internucleon interactions, in 
methods to calculate properties of many-nucleon systems, and in the ability to effectively use grow- 
ing computational power [8]. Because these approaches will be limited in scope for the foreseeable 
future, a natural goal is to develop ab initio DFT for nuclei. In this context, "a6 initio" is taken to 
mean a formalism based directly on a microscopic nuclear Hamiltonian that describes two-nucleon and 
few-body scattering and bound-state observables, in analogy to calculations in quantum chemistry or 
condensed matter physics that start from the Coulomb interaction. This contrasts with many nuclear 
EDF approaches [TT] that fit a functional without relying on an explicit underlying Hamiltonian |12j . 
Efforts to construct bridges between ab initio few-body calculations and the largely empirical nuclear 
EDF's are bringing together diverse theorists and formal techniques using insights from other fields. 
The language and formalism differences are a barrier to progress. We hope to lessen this barrier with 
this review by setting up ab initio DFT as an intermediary. 

There are multiple possible paths to ab initio DFT and the optimal choice for describing nuclei is 
not clear. In confronting the limitations of the most widely used conventional Coulomb DFT implemen- 
tations (such as so-called "generalized gradient approximation" or GGA functionals), condensed matter 
physicists and quantum chemists have made extensive developments toward ab initio Coulomb DFT 
based on wave-function methods. We would like to exploit these advances. This means understanding 
what can be borrowed directly for nuclei and where modifications are needed. At the same time, DFT 
based on effective actions may suggest alternative approximations as well as connections to effective 
field theory (EFT). The goal of this review is to outline various strategies that are being adopted (or 
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may be explored soon), identify common features and challenges, and generally make them more acces- 
sible to the various communities of nuclear physicists attacking these problems. We restrict ourselves 
to a definition of ab initio DFT that is consistent with usage in Coulomb systems (see Section 11.21) 
but which is a subset of the full range of efforts pointing toward non-empirical nuclear EDF (e.g., see 
Refs. [Ml [a [151 [16]). 

The focus on ab initio DFT does not mean we propose abandoning the successes of the empirical 
EDF's, which already achieve an accuracy for known nuclear masses that will be hard to reach directly 
with ab initio functionals. Furthermore, it will only be possible in the near future to make ab initio 
calculations of a limited subset of all nuclei. DFT was originally formulated and is still typically 
described in terms of existence proofs. These proofs imply that it is possible to find a functional (or 
functionals) that depends only explicitly on the density and which is minimized at the ground state 
energy with the ground state density. While these proofs are not constructive, they can be taken to 
justify empirical nuclear EDF approaches. An important prong of the nuclear DFT effort seeks to make 
the EDF's less empirical and therefore more reliable for extrapolation to unmeasured nuclear properties 
by generalizing or constraining the functionals based on ab initio input. This can be done directly 
using constraints from accurate ab initio nuclear structure calculations (e.g., fitting the theoretical 
neutron matter equation of state) but also through insights from ab initio DFT about the form and 
characteristics of the functionals. 

1.2 Basic features/ingredients of DFT 

Our discussion is based on the nuclear many-body problem formulated in terms of a nonrelativistic 
Schrodinger equation for protons and neutrons, with a Hamiltonian of the form 

Hn =f + V = f+V^^ + V^NNN + . . . , (1) 

where T is the kinetic energy and V is the sum of two- and three- and higher-body forces in a decreas- 
ing hierarchy, which is truncated at three-body forces in the most complete present-day calculations. 
(Effects from relativity and other degrees of freedom are absorbed into the potentials either explicitly or 
implicitly.) Such Hamiltonians are derived in low-energy effective theories of quantum chromodynamics 
(QCD) with varying degrees of model dependence. The development of better Hamiltonians, and of 
many-body forces in particular, is an on-going enterprise [T7]. We emphasize that "Hamiltonians" is 
plural because there is not a unique or even preferred form of the short-distance parts of the potentials 
(the longest-ranged part, pion exchange, does have a common, local form in almost all potentials). 
Contrast this with the electronic case, where the long-range Coulomb potential is for many systems the 
entire story. 

For most of our discussion it is irrelevant whether the Hamiltonian being used results from a sys- 
tematic effective field theory (EFT) expansion [17] or a more phenomenological form [l8j, as long as it 
reproduces few-body observables. What is relevant is that the initial Hamiltonian can be transformed 
(e.g., with renormalization group methods) to maintain observables while making it more suitable for 
particular many-body methods (see below). We argue that transformations to soften the potential will 
be critical in making DFT a feasible framework for nuclei; that is, DFT implies an organization of the 
many-body problem that will not work well with all nuclear Hamiltonians. 

We can classify microscopic nuclear structure methods into two broad categories, wave function and 
Green's function methods. In the former, one solves in some approximation the A-hodj Schrodinger 
equation for the A-body wave function "^{xi,- ■■ ,xa), where Xi is shorthand for all of the variables 
of nucleon i (e.g., Xj, spin, isospin). If the operators are known, this allows the calculation of any 
nuclear observable. Methods in this category include Green's functioiJl] and auxiliary field Monte Carlo 

^Despite the name, GFMC is not a Green's function method in the sense it is used here. 
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(GFMC/AFMC) [SI [20], no-core shell model (NCSM) [2l], and coupled cluster (CC) [22l [23]. The 
computational cost of such calculations rises rapidly with A. Nevertheless, most of the recent progress 
in ab initio nuclear structure physics has come from pushing these techniques to higher A UT] [23] . 

Our ab initio DFT discussion will connect to wave-function based formulations that use a single- 
particle basis and can handle non-local interactions (e.g., NCSM and CC but not GFMC). In principle, 
such a formulation solves the problem of finding the ground-state energy Egg of a given (for a 
specified number of nucleons A) by minimizing (^^li^Tvl^) over all normalized anti-symmetric A-particle 
wave functions [1]: 

= min(^l^ivl^) • (2) 

In practice, of course, the Hilbert space (i.e., the basis size) is finite and Egg is found approximately. 
(The calculation is also not variational in many cases, such as CC, but that is not an issue here.) 

An alternative to working with the many-body wave function is density functional theory (DFT) 
[H [H [26], which as the name implies, has fermion densities as the fundamental "variables". We will 
start with DFT as it is typically introduced, citing a theorem of Hohenberg and Kohn (HK) [27] : There 
exists an energy functional Ey[p] of the density p(x)Il labeled by a (static) external potential fext(x) 
such that 

E,[p] =F[p]+ [ rfxt;ext(x)p(x) , (3) 



which is minimized at the ground-state energy Egg with the ground-state density pgs(x). An example 
of Vext is the electrostatic potential from ions in atoms and molecules, as in Eq. f[30l) . The functional F, 
often designated Fhk in the literature, is independent of the external potential Vext and has the same 
form for any A. In this sense it is said to be universal. The HK theorem offers no help in constructing 
F, but is useful in that it gives a license to search for (or guess) approximate energy functionals. This 
would serve as justification for nuclear EDF's except for the disquieting feature that there is no f ext for 
self-bound nuclei, which makes the meaning of Eq. ([3]) unclear. 

A constrained search derivation |28j is a more illuminating alternative to the proof-by-contradiction 
approach to DFT used by Hohenberg and Kohn. We start as before with the minimization in Eq. 
adding an external potential, but now we separate the minimization into two steps [1]: 

1. First minimize over all \E' that yield a given density p(x): 

min(^|#^ + Kxtl^) = min(^|f + + / rfxt;ext(x)p(x) , (4) 

where V is the full internucleon interaction and Kxt = / '^xfext(x)p(x). Define F[p] as the resulting 
contribution of the first term: 

F[p] = (^f°|f + V^|^f") . (5) 

2. Then minimize over p(x): 

E = minE^lp] = min{F[p] + / (ixt;ext(x)p(x)| . (6) 
p P J 

where the external potential fext(x) is held fixed. 
In principle one works at fixed A by introducing a chemical potential p: 

5{F[p] + j rfxt;ext(x)p(x) -p j rfxp(x)} = , (7) 

^As discussed below, we will have multiple densities in practice but considering the fermion density only suffices for 
now. 
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which imphes 

^ + t;ext(x)=/i. (8) 

The chemical potential is adjusted until the density p resulting from solving Eq. ([8]) yields the desired 
particle number A. Or one minimizes only over p that satisfy j rfx p(x) = A. 

As noted by Kutzelnigg [29], the presentation of the HK theorem as an existence proof is often 
accompanied by misleading statements such as "all information about a quantum mechanical ground 
state is contained in its electron density p" or that "the energy is completely expressible in terms of 
the density alone." These claims seem at odds with the observation that while the external potential 
energy is expressible in terms of p (if there is a fext), the kinetic energy is given in terms of the one- 
particle density matrix and interaction energies require two-particle and higher density matrices. How 
do we reconcile this? The key is that the usual wave-function treatment of the many-body problem as 
in Eq. ([2]) has in mind a single, fixed Hamiltonian. In that case, to make a variational calculation of 
the ground-state wave function \E', the energy E must be made stationary to variations in the relevant 
density matrices and not just the density. This corresponds to variations of the normalized A-hodj 
wave function \E'. 

To understand DFT we should consider instead a family of Hamiltonians H[v], each characterized 
by a potential v for which we know the corresponding ground state energy E[v]. We might ask, if we 
know E[v], why not just evaluate at f = fext and avoid more complications? But if we do know E[v], 
then we can construct the functional Legendre transform [29], 



— F[p] = min| J (ixt>(x)p(x) — E[v]^ 



(9) 



where the minimization is over an appropriate domain of v (really the infimum or greatest lower bound) 
rather than just considering a fixed fext [29]. Thus we obtain the dependence of the internal energy 
on the density. This justifies the suggestive notation of Eq. ([3]): If f (x) is set to a constant, it acts 
as a chemical potential and the equation expresses a Legendre transform between two thermodynamic 
potentials. (An expanded version of the thermodynamic analogy is given in Section HI) If the Legendre 
transform is possible (see Ref. [301 [29]), we can also obtain with a second Legendre transformation that 



E[v] = min| J (ixt>(x)p(x) -|- F[p]^ = min|i?t,[p]} , 



(10) 



which is the energy for fixed v expressed as a minimization over a trial set of densities. Thus we 
reproduce Eq. (jH]) and show the origin of the HK expression in Eq. ([3]). 

This perspective shows that DFT and Eq. are really in the spirit of the other major category 
of microscopic nuclear structure methods, namely Green's functions. Instead of the many-body wave 
function, the Green's function approach considers the response of the ground state to adding or removing 
particles [311 [321 [33] . The underlying idea is that knowing the most general response of the ground state 
(or the partition function in the presence of the most general sources) gives a complete specification of the 
many-body problem. Observables such as the ground state energy, densities, single-particle excitations, 
and more can be expressed in terms of Green's functions (with one-body operators needing the single- 
particle Green's function, two-body operators generally needing the two-particle Green's function, and 
so on). The special case considered here starts with the response of the energy to a "source" f(x), 
which is coupled to the density. Instead of non-local sources that individually create particles in one 
place and destroy them elsewhere, here the perturbation by the source is a local shift in the density. 
Because it is a more limited response than the usual Green's functions, the corresponding observables 
probed are also limited, but include the ground-state energy. A natural mathematical framework for 
such responses and other Legendre transforms is the effective action formalism using path integrals. We 
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Figure 1: Kohn-Sham DFT for a fext = harmonic trap. On the left is the interacting system and 
on the right the Kohn-Sham system. The density profile is the same in each. 

give more details in Section H] on how this works. (This perspective also shows that taking f = is not 
a problem in principle; such sources are usually set equal to zero at the end, although in this case there 
are related issues with self-bound systems such as nuclei.) 

In practice, DFT is rarely implemented as a pure functional of the density, such as a generalized 
Thomas-Fermi functional, because no one has succeeded in constructing one that yields the desired 
accuracy (an immediate problem is finding an adequate functional of density only for the kinetic energy). 
Instead, the most successful procedure is to introduce single-particle orbitals that are used in what 
appears to be an auxiliary problem, but which still leads to the minimization of the energy functional 
for the ground-state energy Egg and density pgg. This is called Kohn-Sham (KS) DFT, and is illustrated 
schematically in Fig.[T]for fermions in a harmonic trap. The characteristic feature is that the interacting 
density for A fermions in the external potential Wext is equal (by construction) to the non-interacting 
density in another single-particle potential. This is achieved by orbitals {0j(x)} in the local potential 
f KS ( [p] , x) = t>Ks (x) , which are solutions to 

[-VV2m + t;Ks(x)]0i(x) = eM^) , (11) 

and determine the density by 

A 

p(x) = 5^n#,(x)r = 5^|0,,(x)r, (12) 

i i=l 

where the sum is over the lowest A states with rii = 6{ep — Si) here. When we include pairing, the 
sum is generalized to be over all orbitals with appropriate occupation numbers (see Section [5?T|) . The 
magic Kohn-Sham potential wks([p],x) is in turn determined from 6E^[p]/6p{'x) (see below). Thus the 
Kohn-Sham orbitals depend on the potential, which depends on the density, which depends on the 
orbitals, so we must solve self-consistently (for example, by iterating until convergence). We will return 
later to address the meaning of the KS eigenvalues e^. The ground-state energy is Egs = -Et,^^Jpgs]. 

We will define orbital-based density functional theory (DFT) broadly as any many-body method 
based on a local ("multiplicative") background potential (what we called vks above) used to calculate 
the ground-state energy and density of inhomogeneous systems in the manner just described. That is, 
there will be a single-particle, non-interacting component of the problem that involves solving for orbitals 
with a local (diagonal in coordinate space) potential. We will also require that there are no corrections 
to the density obtained from these occupied orbitals (as in usual Kohn-Sham DFT). It is not obvious at 
this point that this is a necessary feature, because it is not essential to the numerical simplicity or good 
scaling behavior. We will see in later sections how it arises. This characterization of DFT can be realized 
in seemingly very different approaches, such as a particular organization of (possibly resummed) many- 
body perturbation theory (MBPT) and effective actions for composite operators (based on functional 
Legendre transformations). The DFT formalism is often said to be a mean- field approach because of 
the Kohn-Sham potential and this applies to our general definition as well. The point is that it is not a 
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mean-field approximation but an organization that takes a mean-field state as a reference state, which 
if solved completely includes all many-body correlations. (The real issue is how much correlation is 
included in a given approximation to the exact functional.) 



Kohn-Sham Potentials 



Energy 
Functional 



Schrodinger 
-Eqn. Solver 



Orbitals and Occupation Numbers 



Figure 2: Generic self-consistency cycle for Kohn-Sham DFT. The energy functional takes orbital wave 
functions and eigenvalues (with occupation numbers) as inputs. The outputs are the local Kohn-Sham 
potentials from the functional derivative of the energy functional. These could be directly evaluated as 
with a Skyrme or density matrix expansion functional, or solved from OEP equations. 

Implementations of orbital-based DFT will have a self-consistency cycle of the form shown in Fig. [2j 
The code that solves the KS single-particle Schrodinger equations ( "Schrodinger-Eqn. Solver") can be 
generic (if generalized to include pairing) because the same equations are solved for different energy 
functionals (only the potentials change). This part of the calculation is generally the key to the com- 
putational scaling because the cost goes up gently with A and it also means that we can adapt the 
well-developed tools used for Skyrme calculations. The "Energy Functional" will be particular to the 
implementation, ranging from simple function evaluation to solving complicated integral equations. If 
the cost of evaluating this box, which includes calculating the functional derivative defining vks, can be 
kept under control, the computational advantage will hold. (Conversely, in full calculations of orbital- 
based DFT, one has to consider seriously the computational scaling when comparing to alternative 
strategies.) It is often the case that the energy functional is taken to be of the local or semi-localf| form: 

E[pi,T„ ■■■]= j d:>c£{pi{x.),Ti{x.), . . .) , (13) 

where we have allowed different types of densities such as the kinetic energy density and where depen- 
dence on gradients of the densities is also allowed. For example, the Skyrme functional in Eq. fllOll) has 
this structure. We emphasize that this is not a general form, but requires significant approximations to 
be derived from a non-local orbital-based functional (such as by applying the density matrix expansion, 
see Section [5^ . 

We also emphasize that the Kohn-Sham potentials vks are always local no matter how non-local the 
energy density becomes. This is different from some nuclear EDF approaches that feature finite-range 
effective interactions in the form of a Hartree-Fock functional (e.g., Gogny), for which the single-particle 
equations have non-local exchange potentials [TT]. While the locality of Kohn-Sham potentials eases 
the computational burden, it is a constraint that may ultimately prove to be too limiting for nuclear 
energy functionals [T3| [16]. 

We can separate our subsequent discussion of ab initio DFT into two parts: 

"Semi-local" in this context means that the energy density at x depends only on the electron density and orbitals in 
an infinitesimal neighborhood of x [34j . So £ has only a finite number of gradients. 
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1. Given an energy functional, what is the associated Kohn-Sham potential? 

2. How do we construct an ab initio energy functional systematically? 

The second question is more fundamental but the first one is more generic, so we begin with it in 
Section [21 two approaches to the second question are described in Sections [3] and HI We stress that 
Kohn-Sham DFT, with orbitals, is in fact a natural development in each of these approaches, rather than 
simply a "trick" to better approximate the kinetic energy. Further, we find that keeping the densities 
fixed entirely from the Kohn-Sham orbitals either follows as a consequence of the DFT minimization 
conditions or can be used as an imposed condition to derive a functional. In our working definition of 
DFT we allow complete freedom in defining the Kohn-Sham system, which permits more physics to be 
shifted into the potential ("mean field"), making the DFT functional more effective (e.g., a perturbative 
expansion will converge more rapidly). One way of doing this is to allow any local densities paired with 
corresponding sources. 

1.3 Coulomb vs. nuclear DFT 

We will rely heavily on the progress made in ab initio DFT for Coulomb systems, but we should 
always keep in mind the differences between Coulomb and nuclear many-body problems, which will 
introduce substantial challenges. For most Coulomb applications, the Hamiltonian is well known and 
takes a simple, two-body local form (that is, it is diagonal in coordinate representation and has no 
spin dependence). While in principle one could modify the interaction at short distances, e.g., with 
unitary transformations, the original local form is clearly preferred. Because the interaction is to good 
approximation 1/r, it does not make sense to transform it. 

The 1/r potential follows in a straightforward manner from the underlying theory of quantum 
electrodynamics (QED)0 A directly analo gous ab initio calculation of the strong interaction would have 
to start with the quark and gluon interactions of quantum chromodynamics (QCD). But because quarks 
and gluons are not efficient low-energy degrees of freedom because of confinement, low-energy effective 
theories of QCD are used to construct interactions between protons and neutrons. These interactions 
may be systematic (e.g., using EFT) or more phenomenological. But as effective interactions they are 
not unique and transformations may result in forms more amenable to the DFT formulation. 

The weak strength and long range of the Coulomb potential means that the binding energies of atoms 
and molecules are numerically dominated by the Hartree contribution |6J. This dominance would make 
the problem simple except that while the exchange-correlation energy (what is beyond Hartree) is a 
often a small fraction of the total binding energy of atoms, molecules, and solids, it is of the same 
size as the chemical bonding or atomization energy [1]. Thus an accurate DFT functional for this 
contribution (called E^c) is essential and this is the principal challenge of Coulomb DFT, although 
there are complications in certain electron systems (e.g., from pseudopotentials, relativity, etc.). 

It might be supposed that the dominance of Hartree- Fock contributions to the energy, so that cor- 
relations are small corrections treatable in (possibly resummed) perturbation theory, is an important 
reason why Coulomb DFT works so well. In contrast, for typical realistic NN interactions, correla- 
tions are much greater than the Hartree-Fock contribution! This may mean that DFT fails for these 
Hamiltonians. But we can use the possibility of modifying the interactions by renormalization group 
methods [23 EB] (and related methods [H7', fSH]) to change the "perturbativeness" of nuclei. There are 
several sources of nonperturbative physics for typical nucleon-nucleon interactionsjfl 

1. Strong short-range repulsion ("hard-core" for short). 

^If heavier atoms are being considered there will be in practice a more involved potential and possibly the need to 
treat the system fully relativistically. 

^There is also pairing, which we consider separately. 
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2. Iterated tensor interactions (e.g., from pion exchange). 

3. Near zero-energy bound states (e.g., the deuteron and near bound state in the ^Sq channel). 

However, the first two sources depend on the resolution (i.e., the degree of coupling to high-energy 
physics), and the third one is affected by Pauli blocking. Thus we can use the freedom of low-energy 
theories to simplify calculations by lowering the resolution, which softens the potential and makes 
the nonperturbative nuclear physics more perturbative. That is, the convergence rate of perturbative 
expansions such as the Born series for free-space scattering or the in-medium sum of particle-particle 
ladder diagrams improves at lower resolution. This is demonstrated using a quantitative measure of the 
convergence in Ref. [39j . 

In general, transformations to low-resolution interactions that are compatible with the natural reso- 
lution scale of nuclear bound states make the nuclear problem look in some ways more like the Coulomb 
problem. Not entirely, of course, but generally much more perturbative with weaker tensor correlations. 
For example, Hartree-Fock plus second-order many-body perturbation theory may be well converged 
for the bulk energy of the uniform system [39i |40j . Residual issues include how to handle contributions 
from mixing into the ground state of low-lying excitations in finite nuclei, which could require a nonper- 
turbative treatment. In addition, we might expect that such contributions to the energy functional are 
significantly non-local, which may be why empirical nuclear EDF's (which are semi-local) have problems 
incorporating this physics and typically rely on additional procedures to handle these corrections [41j. 
Testing whether orbital-based DFT can accommodate this physics will be an important area for study. 

The conventional DFT for Coulomb systems takes as its starting point precision calculations of 
the uniform electron gas. These are ab initio numerical calculations (e.g., using GFMC) combined 
with well-controlled analytic limits. This solid starting point for a local density approximation (LDA) 
is combined with constrained gradient terms to construct semi-local functionals (GGA) that are not 
microscopic but are also not fit to dataj^ For any DFT, the uniform system is special, because E vs. p 
is a limiting case for the functional that is an observable function. (In general, the functional evaluated 
with a density that is not at the minimum is not an observable, see Section I^TTl ) At present, ab initio 
nuclear calculations of the uniform system with equal numbers of protons and neutrons ("symmetric 
nuclear matter") are much less controlled than the electron gas counterpart. As a result, present-day 
nuclear EDF's are often purely empirical (e.g., most Skyrme interactions): a parametrized form has 
constants fit to properties such as binding energies and radii of a set of nuclei. Indeed, the most 
stringent constraints on the saturation properties of symmetric nuclear matter come from taking the 
uniform limit of the fitted EDF's. Other nuclear EDF's take guidance from the best possible nuclear 
matter calculations but then fine tune to yield quantitatively accurate properties of finite nuclei (e.g., 
see Refs. [HJ |15] and references cited therein). 

The calculation of nuclei also has features without direct parallel in the basic Coulomb systems, 
which complicates the formulation of DFT. Perhaps most apparent is the role of symmetries. Atoms 
and molecules can generally be treated in Born-Oppenheimer approximation, where the slow nuclear 
degrees of freedom provide an external (Coulomb) potential (this is VgxtO for the fast electronic degrees 
of freedom. The basic Coulomb problem for DFT is to find the minimal energy given a configuration 
of nuclei, treated as fixed in space, and the related distribution of electrons (the density or, more 
generally, the spin densities). This external potential means that symmetries of the Hamiltonian such 
as translational and rotational invariance are not realized in the physical ground state. In contrast, 
ordinary nuclei are self-bound and should reflect these symmetries. 

But the nuclear Kohn-Sham potential fKs(x) breaks these symmetries, so that the system being 
calculated is an intrinsic "deformed" state. This is familiar from any mean-field based calculation 

^These are called "non-empirical" (3H [33] because the additional parameters are determined by constraints rather 
than data; quantum chemists also construct empirical Kohn-Sham DFT functionals that are less microscopic. 
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of nuclei [IHl HTj; in general, the organization of the problem about a mean-field inevitably breaks 
symmetries. Handling the restoration of symmetries is well known from a wave-function point of view 
through the use of projection. How we should deal with it for ab initio DFT has only recently been 
considered and may require significant developments [121 SHI SHI- The fact that nuclei are self-bound 
presents not only practical problems but conceptual problems, because there is no external field in the 
case of interest. In the face of symmetry breaking, is the Kohn-Sham approach even well defined? 
This issue has been considered by various authors recently and will be reviewed in Section 15. 2[ This 
is a situation where an alternative formulation, in this case using effective actions, lends a somewhat 
complementary interpretation that can suggest different ideas for approximation. 

Bulk nuclear matter is superfiuid and the treatment of pairing is found to be crucial in the accurate 
reproduction of experimental trends in finite nuclei by empirical energy density functionals. For Skyrme- 
type EDF's, pairing has been accomodated by generalizing from zero-range Hartree-Fock equations 
(equivalent to Hartree) to zero-range Hartree-Fock-Bogoliubov equations [50], with empirical fits of 
the pairing parameters. This requires the inclusion of "anomalous" densities, which is not physics 
that arises for atoms and molecules, and the use of local anomalous potentials leads to divergences. 
Bulgac and collaborators [5T1 [52] have clarified the corresponding renormalization issues and developed 
procedures that could be incorporated into ab initio orbital-based DFT (see Section ISTTl) . An alternative 
path has been followed for Coulomb DFT applied to bulk supercondutivity, where non-local Kohn-Sham 
potentials avoid the new divergences. Similarly, nuclear EDF's with finite-range pairing potentials (e.g., 
Gogny) and recent non-empirical pairing based on low-momentum microscopic interactions in non-local 
form [151 HB] have natural high-momentum cutoffs. While we focus in this review on methods with local 
Kohn-Sham potentials, we stress that the best way to incorporate pairing is an open question. 

All of these features of nuclei impact the energy functional at the same level of accuracy as we 
are trying to achieve, so they cannot be ignored. This further motivates a multi-pronged attack to be 
able to have complementary calculations and to cross-check results, as well as the need to be open to 
alternatives to strict orbital-based DFT. 

1.4 Scope and plan of review 

Our target audience for this review spans several communities of nuclear physicists. Those who do 
wave-function based many-body theory, such as coupled cluster, are likely most familiar with second 
quantization formalism and many-body perturbation theory (MBPT). The practitioners of effective field 
theory, who often come from a particle physics background, are generally more fiuent in the language 
of path integrals and effective actions (although probably not in the form directly analogous to DFT). 
The users of EDF's are experts in the language and techniques of mean-field approximations and broken 
symmetries, dealing with pairing and the like. Rather than concentrate on only one of these formalisms, 
we consider both effective action and many-body perturbation theory perspectives as well as discuss 
how to make empirical approaches less empirical. We believe this will make each more understandable, 
help focus on the key issues, and suggest approximations and generalizations. 

However, it is clear that a complete treatment of these various formalisms would require detailed 
reviews on each topic. Because of space limitations, our discussions will necessarily be rather schematic, 
but we will indicate where to find more details in the literature. Fortunately, there are many fine articles 
targeted at the Coulomb many-body problem that can fill in details. We hope that our treatment here 
will make these articles more accessible to a nuclear physics audience (or, more precisely, audiences). 

Our intention is to provide a guide to possible pathways to improved DFT for nuclei based on the ab 
initio ideas that result in local Kohn-Sham potentials rather than to make an assessment of the current 
status of alternative approaches. Therefore we provide limited details on the successes and problems 
with present-day DFT (or EDF) for nuclei, except for pointers to the literature. We also omit various 
topics, including the many developments in covariant density functional theory, which builds upon the 
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phenomenologically successful relativistic mean-field calculations, time-dependent DFT, the superfiuid 
local density approximation (SLDA), and efforts to analyze and improve EDF's without reference to an 
underlying Hamiltonian. At the end we return to provide brief perspectives on some of these topics. 

The plan of the review is as follows. In Section [2], the formalism for orbital-based DFT is derived 
several ways to help clarify its nature. We also present simplifying approximations that have proven 
notably effective in Coulomb applications. In Section [3], the connection to ab initio wave function 
methods is explored through conventional nuclear MBPT and an improved perturbation approach for 
Coulomb DFT. The critical role of low-momentum potentials to make MBPT for nuclei viable and how 
this can be exploited in a semi-local expansion is also examined. The idea of DFT based on Legendre 
transformations as formulated using effective actions is reviewed in Section HJ with connections to effec- 
tive field theory (EFT) for DFT and conventional Green's function methods. The issues of symmetry 
breaking for self-bound systems, pairing, and single-particle energies are addressed in Section [5l along 
with an overview of efforts to make empirical nuclear EDF's closer to ab initio. Section [6] summarizes 
our perspective on the paths to ab initio DFT and points out some alternative routes. 

2 Orbital-based DFT 

In this section, we introduce the general motivation for orbital-based DFT using the experience of 
Coulomb systems as a guide and derive in several ways how to go from an energy functional to mul- 
tiplicative (local) Kohn-Sham potentials within the optimized effective potential (OEP) method. We 
rely heavily on the reviews by Engel, which appears in Ref. [1], and by Kiimmel and Kronik [53j (see 
also Refs. [SU [6j and Kurth/Pittalis in Ref. [55]), but annotate the presentation with comments on 
differences expected for applications to nuclei. We defer discussion of the construction of appropriate 
energy functionals to later sections. To provide a point of comparison for the OEP method, we start 
with a brief review of the Hartree-Fock (HF) approximation in coordinate representation, which is the 
simplest version of the wave-function microscopic approach. While the exchange-only OEP has close 
similarities to HF, there are important distinctions that persist when OEP includes correlations. 

For simplicity, we present formulas as if the microscopic interactions were always finite-ranged but 
local (i.e., diagonal in coordinate representation), spin-isospin independent, and two-body only. We 
caution the reader that for the nuclear case we will have to deal with non-local forces and (at least) three- 
body forces, all with spin-isospin dependence (see Section [3^ for examples of these generalizations). 

2.1 Hartree-Fock in coordinate representation 

The distinction between the usual wave-function description based on many-body perturbation theory 
(resummed or not) and orbital-based DFT is already evident with Hartree-Fock. As noted earlier, HF 
is a natural starting place for any Coulomb calculation and, with low-momentum interactions, the same 
is now true for nuclei (although Hartree-Fock plus second-order may be a fairer comparison). Further, 
most nuclear EDF's have been viewed at least originally as arising from Hartree-Fock calculations of 
effective interactions (which might be zero-ranged as with Skyrme or finite-ranged as with Gogny) 

Hartree-Fock is the simplest approximate realization of Eq. ([2]), with single Slater determinants of 
orbit als (pi: 

I^hf) -^det{0,(x),2 = (14) 

as the wave functions over which the energy is minimized. The spin-isospin dependence of the nuclear 
interaction is critical to the physics but not to the structure of the equations, where it is a distraction. So 
in order to keep the notation from getting cumbersome, through much of this review we let x represent 
the coordinate variable and, when relevant, also the spin-isospin indices. (Similarly, J (ix includes a 
summation over spin and isospin when relevant.) 
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The Hartree-Fock energy in coordinate representation for a local two-body potential ^(x, y) in the 
presence of an external potential is ^47j 



(vj/HFl^lvi/HF) = /"rfxV0l.V0,+^^ /"rfx/"t/y|0,(x)|V(x,y)|0,(y)|2 

i=l i,j=l 

A A 

-^E /rfx/rfy0l(x)0,(y)V^(x,y)</.](y)0,(x) + 5^ [dyvUy)\Uy)\\ (15) 



i,j=l 1=1 

where the sums are over occupied states. Note that this energy functional of the orbitals is non-local 
in that there are integrals over x and y. The minimization of the energy is achieved by a variation of 
Eq. (|T5l) with respect to the 0,: 



(x) 



— ((M/HF|i^|^HF)-5^^, /rfy|0,.(y)p) =0, (16) 



with the HF eigenvalues Ej introduced as Lagrange multipliers to constrain the orbitals to be normalized. 
There are no other subsidiary conditions, such as imposed below in the OEP. The result is the familiar 
coordinate-space equation for 0j(x): 



A 

VV.(x) + J2 / dyV{^,y)<P*iy){<P,iy)U^) - 0,(x)0,(y)} = £,0,( 

.7=1 



2M 

or, after defining the local Hartree potential and the non-local exchange or Fock potential, 

A 



(17) 



rH(x) 



rfyr(x,y)^|0,(y)P = J dyV{^,y)p{y) 



A 

rF(x) = -r(x,y)5^0*(y)0,(x) = -r(x,y)p(x,y), (19) 



we find the non-local Schrodinger equation: 



2M 



V + rH(x)}0,(x) + j dyrF(x,y)0,(y) = £,0,(x) . (20) 



The equations for the Hartree-Fock orbitals are solved with the same self-consistency cycle as in Fig. |2i 
This is more involved than solving the Kohn-Sham equations flTTl) because of the non-local Fock poten- 
tial, but it is drastically simpler than solving for the full wave function. (A typical numerical solution 
method is to introduce a basis, e.g. the harmonic oscillator basis, which reduces the calculation to a 
straightforward linear algebra problem.) 

The coordinate-space HF equations become significantly more complicated with non-local NN po- 
tentials, such as the soft low-momentum NN potentials, and with three-body potentials. With non-local 
potentials, even the Hartree piece is no longer a multiplicative potential. Indeed, working in coordinate 
representation is not particularly natural in this case, which may raise questions about the appropri- 
ateness of the DFT focus on locality [56]. The treatment of non-local, momentum-space potentials at 
HF is outlined below in Section 13.41 

At the opposite extreme, if the interaction is taken to be zero ranged (contact interactions) so that 
V{'x,y) — i> V{x.)6{'x — y), then the Hartree and Fock terms reduce to the same multiplicative form. 
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This is why "Skyrme-Hartree-Fock" is not only equivalent to a Hartree functional, but has the same 
form as Kohn-Sham functionals that include correlations beyond exchange. This is also why a ladder of 
approximations leading to full orbital-based DFT (as described in Section [H]) naturally involves Skyrme 
functionals on the lower rungs. (However, we emphasize that the microscopic functionals described 
below do not assume zero-range interactions.) 

Finally, we comment on the interpretation of orbital eigenvalues. The Kohn-Sham eigenvalues in 
Eq. (fTTl) appear to be analogs of the HF orbital energies in Eq. (!20ll . The latter are well-defined 
approximations to the separation energies: 



The interpretation of the corresponding KS eigenvalues that appear in the orbital-based DFT treatment 
is less clear, as there is no Koopman's theorem to identify the difference of A + 1 and A body ground- 
state energies with the eigenvalues [53]. However, Gorling [57] has shown that differences of Kohn- 
Sham eigenvalues are well-defined approximations to excitation energies. In addition, Janak's theorem 
holds that the ionization potential equal to the chemical potential is given by the highest occupied KS 
eigenvalue [6] . (This follows in most cases simply from the large-distance fall-off of the physical density 
as dictated by the contribution of the least bound particle, whose wave function falls off according to 
its energy.) The viability of a physical interpretation for the KS eigenvalues is an important topic for 
nuclear DFT and we return to it below. 

2.2 Motivation for orbital-dependent functionals 

The motivation given in the literature for orbital-dependent DFT functionals for Coulomb systems is 
based on some characteristic failures of semi-local functionals as well as the need for greater accuracy in 
some applications. Both of these points are relevant to nuclear DFT, where EDF's of the Skyrme type 
are also semi-local and where greater accuracy is sought globally and greater reliability is sought for 
extrapolations. Before considering the motivation in more detail, we first review some of the standard 
Coulomb DFT formalism, indicating where our treatment for nuclei will differ. Following Engel [1], we 
restrict our discussion to non-relativistic, time-independent, spin-saturated systems; this is for simplicity 
and is not a limitation of the formalism. 

The total energy functional is generally decomposed for Coulomb systems as (except that n is 
typically used for density instead of p; other details of the notation also vary in the literature) 



where is the KS kinetic energy, E^^t is the external potential energy, Ey{ is the Hartree energy, 
and E^c ("xc" stands for "exchange-correlation") is defined to be everything that is left over; i.e., it is 
implicitly defined by specifying the other pieces. (Note: We have omitted a piece describing the energy 
of the ions themselves.) For Coulomb systems with just a 1/r potential there are explicit expressions 
in terms of KS orbitals for all but i?xc, namely 




for a > A (particle) 
for a < A (hole) . 



(21) 
(22) 



^tot[p] = T,[p] + Eext[p] + E^\p] + E^,[p\ 



(23) 




(24) 




(25) 




where V"c(x, y 



x-y 



(26) 



13 



As before (and throughout this review), these orbitals satisfy 



•^ + -Ks(x) 



(x) , (27) 



where the KS potential fKs(x) = ^(-E'ext + -^h + -Sxc)/'^p(x)Jll The signature features for ab initio DFT 
(in our broad definition) is that f ks (x) appears multiphcatively in the KS equation and that the density 
is given completely by summing up the occupied (defined as the energetically lowest here) KS states: 

p(x) = ^n,|0,(x)|2. (28) 

i 

At finite temperature or when pairing is introduced, the sum will be extended to all orbitals with 
appropriate occupation numbers rii (see Section ISTTl) . 

Given Eq. ( I23i) . the KS potential fKs(x) is the sum of three pieces, 

vks{^) = ^^cxt(x) + vii{x) + Wxc(x) , (29) 

where 

-Efe. (30) 



WexttXl 



a=l 



5p(x) J |x - x' 



^h(x) = -YZTZTT = ^ / T:: TTi ' (31) 



and 



These formulas for fext(x) and fH(x) are particular to the Coulomb problem, but the structure of these 
terms is more general, as is the expression of fxc(x) as a functional derivative. 
Some comments about these formulas: 

• While real nuclei usually do not have external potentials, it can be useful to theoretically put 
the nucleus in a trap for comparisons of empirical functionals to ab initio calculations. However, 
the external potential should be viewed more generally as a source that is varied and then set 
to zero at the end, such as those used in field theory with path integrals. Therefore we can 
add more general sources coupled to local densities, such as the kinetic energy, spin-orbit, and 
pairing densities found in Skyrme EDF theory and expect to have better energy functionals. In 
this case, we have a series of Kohn-Sham potentials, each equal to a functional derivative of the 
corresponding density. We will see how this works explicitly for the kinetic energy and anomalous 
(pairing) densities in Sections 14.21 and [STTl respectively. 



If a functional in the semi-local form of Eq. f[T3|) is constructed (e.g., from a local density ap- 
proximation plus gradient corrections), then v^c can be evaluated trivially from Eq. fl32l) and the 
self-consistency cycle of Fig. |5] can be carried out directly. This is the form of Kohn-Sham DFT 
that is most widely used (e.g., GGA and its variations for the Coulomb problem and Skyrme, 
SLDA, etc. for the nuclear problem). 

The success of Kohn-Sham compared to Thomas- Fermi (for which the kinetic energy is a functional 
of the density) followed from the introduction of orbitals to treat the kinetic energy, with manifest 
improvements such as the reproduction of oscillations from shell structure [26]. Thus, in this sense 



^The KS potential, density, etc. are often denoted Vg, Ps, and so on in the Coulomb DFT literature. 
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orbital-based functionals are already inevitable and making the potential energy orbital dependent 
is not a major step conceptually [4J. Note that the kinetic energy contribution here is just a 
definition and, while a good approximation, it is not the many-body kinetic energy (which would 
be expressed in terms of the full single-particle Green's function or the one-body density matrix). 
The exchange-correlation part has to absorb the difference between these kinetic energies. 

• For Coulomb systems the Hartree contribution is dominant and the Hartree energy functional 
takes a special form that depends explicitly on the density only. These are the reasons for singling 
it out. For nuclear systems the Hartree piece is neither the dominant contribution nor simple in 
general with low-momentum potentials, which are non-local and require integrals over quantities 
that do not reduce to the density (e.g., the one-particle density matrix). So singling it out is not 
generally helpful. It is advantageous in some approximations such as the density matrix expansion 
(DME) to isolate the long-range part of the potential, which is local, and to treat this part of the 
Hartree potential explicitly. 

• The exchange part can also be singled out for Coulomb because of the possibility of a precise 
formula in terms of the KS orbitals and the Coulomb potential. But isolating particular parts of 
the functional is not necessary, because 

t;H(x) + t;.c(x) = t;Hxc(x) = -^{En[p] + E^,{p\] ^ , (33) 

op(xj op(xj 

so all the interaction terms can simply be combined (and we use the "Hxc" subscript to indicate 
when this is done). This is the form we generally use for the application to nuclei. 

As is well documented, the conventional Kohn-Sham DFT with semi-local functionals has been 
extremely successful for Coulomb systems. However, there are failures and these help motivate the 
development of orbital-dependent functionals. It is not yet clear to what extent these failures have 
analogs for nuclear-based system, but we summarize the list of arguments (taken largely from Engel's 
review in Ref. [1]; look there for specific references) along with some speculations about nuclei: 

• Heavy Elements. For heavy constituents (e.g., gold), the local density approximation (LDA) 
tends to work better than the generalized gradient approximation (GGA) and this is not at- 
tributable to relativistic effects. The suggestion is that GGA has trouble with higher angular 
momentum {d and /). If this is a generic problem it would certainly be relevant for nuclei. The 
fact that GGA is not a systematic improvement over LDA is also motivation for ah initio DFT 
— to develop a hierarchy of approximations that do systematically improve, as for the coupled 
cluster method 



Negative Ions. The fall-off of the Kohn-Sham potential does not have the 1/r asymptotic form 
needed for negative ions and Rydberg states. The physical picture is that if one electron is far 
away from the others, it should see the net charge of the remaining system of — 1 electrons and 
protons. But because f h still always has the Coulomb repulsion of the far electron, it has to be 
removed by fx, but this does not happen with LDA or GGA functionals. Engel emphasizes that 
the exchange functional needs to be rather non-local to cancel the self-interaction in the Hartree 
potential [1] . For the nuclear case, the impact of the self- interaction problem was considered long 
ago in Ref. [60] but only recently reconsidered along with the additional problem of self-pairing 
in Refs. 



Dispersion Forces. Dispersion forces are a type of van der Waals force. The problem here is the 
locality of the exchange-correlation functional. If two atoms are so separated such that there is no 
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overlap in the densities, then the density is the sum of the two atomic densities. But we expect 
virtual dipole excitations leading to molecular bonding (this is called the London dispersion force). 
This does not work for LDA because the binding energy from the correlation functional is 

E, = E^^"^ [p^ + ps]- E^^"^ [p^] - [ps], (34) 

which vanishes because only regions with non-zero density contribute to the correlation energy 
(so the first term on the right side is the sum of the other two terms). The same result holds for 
GGA because the density only in the near vicinity of x contributes to the energy density at x. So 
we need non-locality for virtual excitations. Analogous issues for nuclei may arise from coupling 
to low-lying vibrations. 

• Strongly Correlated Systems. Here there are failures for certain systems, such as some 3d 
transition metal monoxides, which the LDA and GGA either predict are metallic when they are 
actually Mott insulators or else greatly underestimate the band gap. Indications are that the 
incorrect treatment of self-interaction correction is the problem (although this is not proven). 

The self-interaction problem can be illustrated by considering the "exact exchange" functional 
of DFT, which is defined to be the Fock term as in Eq. (fT5|) written with KS orbitals. So for a local 
potential V^(x, x'), this is simply 

= Xl^'^^' /^^ / d^'(f)i{^)M^)Vi^,x')(f)l{x')M^') ■ (35) 
kl J J 

This is not the usual HP exchange contribution, because while it agrees in form with Eq. (fTSll . it does 
not have orbitals that satisfy the non-local HP equations. Just like the difference in the kinetic part, 
the difference between HP and KS exchange is absorbed into the correlation functional by construction. 
The functional E^ is a density functional in that the orbitals are uniquely determined by the density, 
but it is imphcit in the density dependence. The more general -Ehxc[p] will also depend on the KS 
eigenvalues. 

The form of E^ ensures exact cancellation of the self-interaction energy in E-^ jl] . Suppose we split 
E^ into two pieces according to whether or not k = I in the double sum. Then 




Note that the integrand in the second term does not reduce to the product of densities because there is 
only one k sum. However, this second term does cancel the corresponding part of the Hartree functional 
if we rewrite the latter in a similar form: 




This cancellation is familiar from ordinary Hartree-Pock. But if E^ is expanded in semi-local form as in 
the LDA or GGa for Coulomb systems or the DME for the nuclear case (see Sect. 13.41) . this cancellation 
is lost. 
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Kiimmel and Kronik [53] emphasize the formal deficiencies that can lead to qualitative failures of 
the LDA and GGA predictions. These include not only the non-cancellation of self-interaction but the 
lack of a derivative discontinuity in i^xc- The latter issue starts with the DFT definition of a chemical 
potential /i: 

(5p X 



which is position independent when evaluated at the ground-state density. Perdew et al. argued [6l] 
that this chemical potential must have a discontinuity: if the integer number of electrons is approached 
from below its absolute value should be the ionization potential while if from above it should equal the 
electron affinity. The derivative discontinuities at integer particle numbers in general comes from both 
the noninteracting kinetic energy and E^^.- But the LDA and GGA -Exc's are continuous in the density 
and its gradient, so there is no particle number discontinuity. Ref. [S3] has more details on why this is 
an important issue for Coulomb DFT. We know of no investigation, however, into its impact on nuclear 
EDF's. 

In the end, a basic issue is that any full DFT energy functional must have non-localities and it may 
be problematic for nuclear structure to expand in a semi- local functional [T2t HHl HQ] . The only way we 
know to test this is by comparing such functionals derived from microscopic interactions (e.g., through 
the density matrix expansion, see Section 13. 4p to a full orbital-based functional. A program to carry 
out such comparisons was recently initiated as part of the UNEDF project [H [9]. 



2.3 Derivation of the optimized effective potential 



The fundamental problem in extracting the Kohn-Sham potential is to calculate the functional derivative 
with respect to the density as in Eq. ([3SD- (More generally, we will need to take functional derivatives 
with respect to other densities, such as the kinetic energy density and the anomalous density for pairing 
for nuclear applications.) Starting from an ah initio energy functional, one way to proceed is an 
expansion such as the density matrix expansion (DME), which results in functionals with the densities 
explicit (see Section [3^ . In the absence of such approximations, the density is naturally explicit only 
in the Hartree energy functional (and in this case only for local interactions). Therefore, we need an 
implicit calculation of the derivatives [1]. 

The most direct procedure is to use the chain rule [1]. As reviewed in later sections, the energy 
functional can be built from the KS orbitals and eigenvalues (e.g., with Kohn-Sham MBPT as in 
Section [3]), so we express the density derivative in terms of those. As an intermediate derivative we 
vary the KS potential vks- 



5p(x) 



dx.' 



c/x' 



5p(x) 



5p(x 

idx" 



^l(x") 



SE^ 



+ c.c. 



+ 



dE^ 



(39) 



Note that in all of the expressions given here the sum over orbitals k is not restricted to occupied states 
unless Uk explicitly appears. (We assume filled shells here and refer the reader to the cited literature for 
the case of unfilled shells.) Now we need to find the functional derivatives introduced in this expression j§ 
For the exchange functional defined in Eq. fl35|) with a local potential, we obtain 



SE, 



64i 



Y^niM^') /"cix0Kx)0fc(x)y(x,x') 

I 



(40) 



^In general one also needs to vary the occupation numbers, although this is not needed if there is fixed particle 
number i62i. 
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and dE^/dsk = 0. The functional derivatives 50|./5fKS and Sek/Sv^s follow from simple first-order (non- 
degenerate) perturbation theory with 6vks as the perturbation. That is, we treat the KS equation for 
the orbitals just like in conventional quantum mechanics. So we start with the change in the eigenvalue: 

SsklSvKs] = J rfx' 0l(x') SvKsi^') M^') , (41) 
which directly implies the functional derivative 

-0l(x)0,(x). (42) 



Similarly, for the wave function, 

6M^) = E f 0f (x")'5^^Ks(x")0.(x") , (43) 

777 — J 



and 



^4(x) = E / 0l(x")5t;Ks(x")0Kx")^^ , (44) 
777 J £k — £l 



l^k 

which imply the functional derivatives 

^M^) V^0KX)0[(XO , / A ^ / , / A /.P-N 

T = 0fc(x) = -C;fc(x,x)0fc(x) , (45) 

dt;Ks(x') £k-£i 

and 

i^ = 4M5:^ipM_4MG.(x'.x), (46) 

where (note the sign change in the denominator, which can differ in the literature) 

G,(x,x').y::«^l^. (47) 

5? 

Note that to construct the Green's function for every state requires explicit solutions for all the orbitals. 
Below we discuss approximations that avoid the complication of unoccupied orbitals. 

To complete the chain rule we need to evaluate 6vks/^P- To do so, we focus first on the inverse 
function, which is the static response function of the KS system, 

= x.(x,x') = -5^n,4(x)Gfc(x,x')0fc(x') + C.C. . (48) 

This follows directly by applying our functional derivative formulas to the expression for p(x) in terms 
of orbitals. Note that only the terms with / unoccupied in will contribute to Xs-, because those with 
/ occupied will cancel (e.g., pick an / and a k and note that the complex conjugate term is the same 
but with the opposite sign energy denominator). 

By multiplying 6Exc/Sp{x.) by Xs and integrating over x, the OEP or 0PM (optimized potential 
method) integral equation is obtainedl^ 



j (ix'xs(x,x')t;xc(x') = Axc(x) 



(49) 



^The names OEP and 0PM are used in the literature by different authors to describe either this equation or more 
generally the orbital-based approach. 
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where 



A, 



dx.' 



^I(x)G,(x,x') 



+ c.c. 



(50) 



This is a Fredholm integral equation of the first kind. Note that because Eqs. (H9l) and (1501) are hnear 
in £^xc5 we can consider separate equations for different pieces of E'xc- It is common to introduce the 
"orbital shifts" [53| 



t/ 



rfx'0|.(x') [Mxc,fc(x') - t;xc(x')]G'fc(x',x) 



where 



SO that we can write the 0PM integral equation compactly as 



()fe(x) + c.c. 



. 



(51) 
(52) 

(53) 



In carrying out the self-consistency loop (Fig. [2]), the solution for the orbitals when given a Kohn-Sham 
potential proceeds the same as always. Given the new orbitals and eigenvalues, we first find (^^(x, x') 
and u^ck, and then the new KS potential by solving the 0PM equation, and the loop repeats. As 
discussed in Section \2.5\ this is numerically difficult and comparatively inefficient, but there are also 
good approximations that simplify the solution significantly. 



2.4 OEP from total energy minimization or density invariance 

To get more insight into the physics of the OEP, we turn to the original derivation of the 0PM integral 
equation, which is based on the minimization of the energy functional with respect to the density [1]. 
Without explicit dependence on the density we do not know how to do that minimization, but the 
Hohenberg-Kohn theorem that tells us that vks and p are directly related. In particular, this implies 
that we can replace the minimization with respect to p by a minimization with respect to f ks (for fixed 
particle number), 

6Etot [4>k,ek] 







(54) 



Now we apply the chain rule as before: 



6E, 



tot I 



• £k] 



6E, 



tot 



+ c.c. 



+ 



Ssk dE^ 
5t;Ks(x') dsk 



(55) 



The functional derivatives of -Etot with respect to the orbitals and eigenvalues are given by 



5E, 



tot 



dE, 



tot 



dei 



dE- 



2M 



+ t;ext(x) + t;H(x) 



>fe X 



dei 



SE^ 



(56) 
(57) 



with the remaining derivatives given by previous expressions. The first term on the right side of Eq. ([56 
can be rewritten because 0a:(x) satisfies the KS orbital equation. 
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Plugging everything into Eq. fl55l) gives the minimization condition: 



J2 [ rfx'(4(x)G,(x,x') \n,M^){v.c - e,) + + c.c.) + \M^)\'^ = . (59) 

Finally, the 0PM integral equation is recovered after identifying Xs(x, x') and Axc(x) and using 

j rfx0l(x)G,(x,x') = j rfx'Gfc(x,x')0l(x') = . (60) 

(Note that Eq. flBUl) means there will be complications with inverting Xs-) 

If we just include the direct and exchange terms in the functional then i?tot looks just like an 
HF functional. The key difference is that the HF approach corresponds to a free (unconstrained) 
minimization of the total energy functional with respect to the (pk and Ek- But the minimization here is 
not free; rather, the 0^ and Sk have to satisfy the KS equations with a multiplicative potential. This is a 
subsidiary condition to the minimization of -Etot, which is implemented by the 0PM equation. Because 
this means the variational calculation is over a more limited set of states, the OEP applied to exchange 
only will always give a higher energy than Hartree-Fock [3]. 

The 0PM equation can be derived yet another way, which is based on the point-by-point equivalence 
of the KS and interacting densities [631 EH [65] : 

p,(x)-p(x)=0. (61) 

This is perhaps the least obvious feature of Kohn-Sham DFT. It might appear in the context of a 
perturbative expansion (see Section [3]) to be simply a choice of the Kohn-Sham potential that makes 
higher order corrections to the density vanish. However, in fact it implies the energy minimization that is 
common to all DFT applications, that is, that the KS exchange-correlation potential is the variationally 
best local approximation to the exchange-correlation energy [M]. The basic demonstration starts with 
expressing Eq. fl6Tl) in terms of traces over the KS and full single-particle Green's functions evaluated 
at equal times and the same coordinate arguments: 

- iTr{G,(xt, xt+) - G(xt, xt+)} = , (62) 

where the Green's functions are defined as usual as time-ordered products of the field operators in the 
KS or fully interacting ground states: 

G,(xt,x't') = -z(<l'Ks|T^^(xt)V^^(xr)|<l>Ks) (63) 
G(xt,xY) = -z(^o|TV'(xt)^t(x't')|*o) • (64) 

The specification t"*" serves to order the field operators as ip'^ip. The full Green's function is related to 
the KS Green's function by a Dyson equation [using a four- vector notation x = (x, if:)]: 

G{x,x') = Gsix,x') + J dydy'Gsix,y) Shxc(i/,Z/0 - - y')^Hxc(y) G{y',x') , (65) 

which follows from the separate Dyson equations for Gg and G by forming G~^ and Gj^ and subtracting, 
solving for G~^, and then taking the inverse [4j. Note that the Hartree parts of the irreducible self- 
energy T,iixc{y,y') and the Kohn-Sham self-energy 6{y — y')^Hxc(y) cancel for local potentials, leaving 
the difference of Exc and v^c- This is known as the Sham-Schliiter equation and is a nonlinear equation 
for Vxc, which can be shown to be equivalent to the 0PM equation. We leave the details of showing 
this equivalence to the references, but give a schematic demonstration in Section HI 
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2.5 Approximations 



While the preceding demonstrations illustrate formal constructions that suffice to carry out orbital- 
based DFT (given appropriate energy functionals), solving the equations can be numerically difficult 
and there are also efficiency issues. (For recent work on numerically stable methods to solve the OEP 
equations using Gaussian basis sets, see Refs. [59l[66].) In Coulomb applications, 0PM calculations are 
found to take one or two orders of magnitude longer than corresponding GGA calculations. The source 
of this inefficiency is the need for the Kohn-Sham Green's function, which requires knowledge of all the 
unoccupied as well as occupied orbitals. An approximation (and subsequent variations) by Krieger, Li 
and lafrate (KLI) |67j based on a closure approximation avoids this problem and seems in practice (for 
Coulomb cases) to be quite accurate, at least for the exchange part of the functional. 

The key is to replace the energy denominator in the Green's function by an averaged difference: 
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(66) 



Upon substituting into the 0PM integral equation, one obtains 
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(67) 
(68) 



Consistent with the closure approximation, the derivative dE^ddSk is neglected, which finally yields 



c.c. 
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KLI 



(69) 



Although appears on both sides, one can either iterate to self-consistency starting from (for exam- 



ple) an LDA approximation to Av^^^ or recast these as a set of linear equations allowing Aw 



KLI 



to be 



determined without knowing 

Results from representative calculations comparing 0PM and KLI to LDA, GGA, and HF using 
exchange-only functionals have been studied systematically and are summarized in Refs. [1] and |53j . 
Note that in such comparisons the Hartree-Fock (HF) result serves as the "exact" answer. In most 
cases, exchange-only 0PM is found to be a very good approximation to HF, e.g., agreement at one part 
in 10~^ for the ground-state energies of heavy closed-subshell atoms. Because exchange-only 0PM is a 
more restricted variational minimization than Hartree-Fock, the HF results must always be more bound 
than 0PM, but the small difference implies that the greater variational freedom for HF has little effect 
in practice. Furthermore, the KLI energies are very good approximations to the full 0PM results; for 
this same example the KLI-OPM difference is systematically about 1/3 that of HF-OPM. The KLI 
energies are above 0PM in all CclSGS, clS required because the full 0PM energy is a minimum for local 
potentials. The deviation for even the largest atom is still very small. The level of agreement can be 
calibrated by comparison to LDA and GGA results; the latter are much better than the former, but 
sometimes an order of magnitude worse than KLI and without a systematic sign. 

Other approximations related to KLI but with improvements have been proposed. The Common 
Energy Denominator Approximation (CEDA) [SS] and the Localized Hartree-Fock (LHF) approxima- 
tion [69] are the same as KLI except that only the energy differences for occupied-unoccupied pairs 
of states are approximated by an average, while the exact differences are kept for occupied-occupied 
pairs. CEDA is invariant under unitary transformations of the occupied orbitals, which is a plus, but 
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in practice the results are similar to KLI ^3] . Yet another approximation is the effective local potential 
proposed by Staroverov et al. [7D], which is efficient to implement numerically. 

The good results that have been found using the KLI and similar approximations (for total energies) 
are encouraging because the corresponding potentials are far easier to calculate. KLI has been shown to 
be a sort of mean- field approximation to the full OEP, which accounts for its success [53] . However, more 
work is needed on applications beyond exchange-only functionals. There are no calculations yet that 
compare these approximations (and more extreme approximations such as the density matrix expansion, 
see Section 13.41) for nuclear systems, so drawing conclusions on the prospects would be premature. 

3 DFT and ab initio wave function methods 

The development of accurate ab initio functionals for orbital-based DFT that go beyond exchange-only 
(what is called the "correlation functional") for Coulomb systems is the subject of much recent activity. 
Progress is being made although it is still an open question whether the most ambitious accuracy goals 
can be met (e.g., chemical accuracy). The extensive review by Kiimmel and Kronik [53j provides a 
good idea of the state of the art in 2008 but we note that there are subsequent and ongoing advances. 

For the nuclear problem, the development of analogous ab initio functionals suitable for orbital- 
based DFT is still in its infancy. Thus our "review" in this section will mainly look toward the future, 
considering first how to formulate an exchange-correlation functional in MBPT, neglecting some serious 
issues such as symmetry breaking. We then examine a specific implementation by Bartlett and collabo- 
rators for quantum chemistry that is based on requiring density invariance and which highlights how the 
convergence of the perturbation expansion for the functional can be improved by shifting more physics 
into the Kohn-Sham potential. These developments are based on applying MBPT at low order, so we 
next discuss how this is meaningful for nuclear interactions if they are transformed to low-momentum 
potentials. Finally, we review recent and ongoing efforts to apply the density matrix expansion to make 
ab initio energy calculations by expanding functionals about the uniform system (asymmetric nuclear 
matter). This also provides examples of how to deal with momentum space potentials, non-localities, 
and three-body forces, which will not be found in the Coulomb literature. 

3.1 Goldstone many-body perturbation theory 

A direct method to construct an orbital-based energy functional is to consider the Kohn-Sham potential 
vks(x) as defining the single-particle potential (typically called f/(x) in the nuclear physics literature) 
that is used in Goldstone many-body perturbation theory (MBPT). This will give us a diagrammatic 
expansion for the energy as a functional of U that will depend on KS orbitals and eigenvalues. There are 
various alternative ways to formulate MBPT. For example, in the Coulomb DFT literature a formalism 
for perturbation theory based on coupling constant integration [H Hj is typically used. Here we use the 
Goldstone formalism that historically led from "hard core" NN potentials to the hole-line expansion [Til 
[72| [73] . but the end result is basically the same. Note that "perturbation theory" does not exclude 
infinite summations of diagrams in our functional, although we will argue in Section 13.31 that second- 
order perturbation theory is a good starting point for low-momentum nuclear interactions. 

Our schematic discussion is compatible with the detailed, pedagogical expositions of the Goldstone- 
diagram expansion given in Refs. [TH [TH [221 [Z3] • We start with a division of the second-quantized, 
time-independent Hamiltonian: 

H = Ho + Hi, (70) 

where the separation is our choice. One might imagine taking Hq to be the kinetic energy plus the fixed 
external potential and Hi to be the interaction potential. Such a division will be invoked in Section 15.31 
when we discuss the connection to the perturbation expansion using Green's functions. However, for 
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Figure 3: Goldstone diagrams at first and second order in Hi with a single-particle potential —U (dash 
lines) and a two-body potential Vnn only (wiggle lines). 

our applications, we will require Hq to be the sum of the kinetic energ}0 and a single-particle potential 
(including any Kxt), Hq = T + U. Then Hi corrects for the new part of U , Hi = V — U + Kxt- With 
this restricted form, the many-body ground- and excited-states of Hq will be simple to construct. The 
ground state of Hq, designated |$o), will be our reference state and our Fock space will be built from 
particles and holes with respect to it. More precisely, the complete set of orbitals defined by the single- 
particle potential establishes the single-particle basis for the Fock space, and second-quantized field 
operators are expanded in this basis. So |<I>o) = n£i<eF'^I|0) a Slater determinant of orbitals where 
ai(aj) is the destruction (creation) operator for the single-particle state i and ai\0) = for all i and 
is the energy of the last filled level. The idea will be to identify U with the Kohn-Sham potential [i.e., 
U = Vks = /'^xfKs(x)p(x)], which is to be determined self-consistently in accordance with Section O 
So the orbitals 0, are then KS orbitals. 

Goldstone's theorem expresses the energy of the ground state as 

- °° / 1 - \" 

E = Eo + {%\HiJ2 [ l'^0)connected , (71) 

where E^ is the ground state energy of Hq\ that is, //qI'^'o) = -^'ol'^o)- (Note: These operators are in the 
Schrodinger picture.) The general idea of an Vb^ order contribution (see Fig. H]) is that we start with 
|$o) and apply Hi^ which creates particles and holes. For a two-body potential, this is two particles 
and two holes (called "doubles" in quantum chemistry parlance), but in our case Hi in addition has 
one-particle-one-hole excitations ( "singles" ) as well as few-particle-few-hole excitations from few-body 
interactions ("triples" and beyond). The factor [Eq — Hq)~^ propagates the state and then the next Hi 
hits. The proviso "connected" means that we do not have |$o) as an intermediate state; we only get 
back to it at the end after the final Hi. 

Goldstone diagrams can be used to organize the contributions from Eq. (ITTll . A set of rules is given 
in Ref. [71] for the general case of interest with nonzero U . The diagrams that contribute in a finite 
system at first and second order in Hi are shown in Fig. [3] [75] (those in the second row vanish in 
uniform matter by momentum conservation). A sample diagram contributing to Eq. fl7Tl) is shown in 
Fig. m Lines with upward arrows are particles and those with downward arrows are holes; they will 
carry labels for the KS orbitals. By inserting complete sets of states, any diagram is reduced to products 
of matrix elements of Hi and energy denominators that are sums and differences of KS eigenvalues. For 
example, a general expression for the energy shift E — Eq with f/ = is found to be (see Ref. [32] for 



For finite nuclei, we would subtract the center-of-mass kinetic energy Tcm from T. 
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Figure 4: A Goldstone diagram contributing to Eq. fl7ip at fourth order in Hi with a two-body potential 
only (wiggle lines) pi] . 

more details): 

connected ^ — ^ — 

where the matrix elements are antisymmetrized Hugenholtz matrix elements (and only the two-body 
interaction is included; the generalization to three-body and beyond is similar only with few-body 
matrix elements such as {ij k\V^T<f^\lmn}) . The energy denominator comes in between each successive 
interaction and includes all particle lines ( "A" ) and hole lines ( "a" ) cut by the dotted lines as in Fig. HI 
The number of hole lines is Uh and the number of equivalent pairs is Ug, while ul is the number of 
closed loops that the diagram with only direct matrix elements at each vertex would have. 

A perturbative calculation applying Goldstone's theorem (which is just a restatement to all orders of 
time-independent perturbation theory for the ground-state energ}0) is organized in the number of times 
Hi is applied. Experience using nuclear potentials with strong short-range repulsion reveals that such 
calculations do not converge. This led to the development of the Brueckner-Bethe-Goldstone approach 
and the hole-line expansion. The latter is a prescription for infinite resummations of diagrams (e.g., 
summing the potential first into a G-matrix and then including all diagrams with a given number of 
independent hole lines |7T]) along with constraints on the background potential U in order to cancel 
certain diagrams. These constraints are incompatible with the choice U = Vks- However, with low- 
momentum interactions we are (apparently) free to choose U in this way and even a simple perturbative 
truncation (at sufficient density) is a reasonable approximation. 

The matrix elements of the low-momentum potentials are typically given in momentum represen- 
tation or, more appropriately for the present discussion, in a harmonic oscillator basis. Then the KS 
orbitals are expanded in the same basis so that matrix elements such as {^jIH^nI^O^ which is in the 
orbital basis, can be evaluated. (Note: the necessary angular momentum recoupling is conventional; see 
for example Ref. [3H].) To calculate fKs(x) according to the OEP prescription, we need to compute the 



^^A perturbative expansion with the same Hq (often called Hg in the DFT literature) and Hi can be derived [1] using 
a coupling constant integration to adiabatically switch on Hi as in Refs. [5T1 132j . Thus H = Hq + XHi and 

E-Eo = Ei^ f dX{^oW\Hi\^oW , (73) 
which is developed using the interaction picture time-evolution operator. Then E^c is given by 

E^c^Ei+ /dxp(x)?;xc. (74) 
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functional derivatives of -Ehxc [see Eq. fl33l) ] with respect to the orbitals and eigenvalues. The eigenvalues 
appear explicitly in the Goldstone expansion and the derivatives of the potential matrix elements can 
be carried out using the basis expansion without having to evaluate the matrix elements in coordinate 
space. For example, the two-body contribution from the Hartree-Fock terms is 



6 ^ 



(x) 1 2 

j,k=l ^ j,k=l 



I Yl lijk\V^^\jk) - {jk\V^^\kj)]\ = J2 M^)[{jk\V^^\jt) - {jk\V^Nm . (75) 



If the matrix element integrals are written out in coordinate space for a local Vnn(x, x') and completeness 
of the orbitals used to remove the sum over k, the corresponding part of Eq. (ITTl) is reproduced. 
So we take [76] 

Ho = H,=f+ /dx p(x)i;ks(x) , (76) 



where vks = Vext + + v^c = Vext + ^Hxc- The KS energy and density are 

Eq = Es=Ts+ j rfxp(x)wKs(x) = ^UiEi , (77) 

p(x) = ($o|p(x)|$o) =5Zri#,(x)|2 = pKs(x) , (78) 

i 

which is the exact density according to the conditions (to be) imposed on wks- (Note: p(x) = ^'''(x)'^(x) 
where the field operators are V'^(x) = ^i4>*i{^)(^\ and ip{yi) = Xli */'«(x)'^l-) Then we have 

H, = V- y"rfxp(x)t;H.c(x) . (79) 

At this point, we can carry out the Goldstone expansion to any order or with whatever infinite sum- 
mation we wish. 

We note that because this expansion depends on fKs(x), it is a highly nonlinear functional equation, 
as £^Hxc depends on its own functional derivative. This is not in conflict with DFT because fHxc is a 
density functional, so 



Eyixc = Eu + E^c = Ei + j dx [i;h(x) + t;xc(x)]p(x) 



^0) 



is an implicit functional. Furthermore, it is the same nonlinearity as in the 0PM equation. A perturba- 
tive approach to solving for the functional derivative fHxc(x) = 5-E'hxc[p]/'5p(x) leads to a linearization 
of the 0PM equation. This will correspond to the inversion method discussed in Section HI The details 
of such a perturbative expansion and the nonperturbative resummation of ring diagrams (the random 
phase approximation or RPA) is described by Engel (with references to the original literature) [1]. 

According to Engel the results for KS perturbation theory are not satisfactory. However, more 
recent developments have shown promise by improving the perturbation theory. We have emphasized 
that with low-momentum potentials, U can be chosen freely to allow the DFT constraints to hold. But 
there is still freedom in defining it, which can be used to significantly improve perturbative convergence 
in the Coulomb problem. We consider this in the next section. 

3.2 Improved perturbation theory 

In this section, we summarize the basic features of the approach to ab initio DFT under development 
by Bartlett and collaborators [581 El] . In Ref . [59] they explain their strategy, which is closely in line 
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with the highly successful coupled cluster approach to electronic systems. A key part of the formulation 
as well as its numerical execution is to use a finite basis set (which is almost always chosen to be 
gaussians for finite systems in quantum chemistry so that all integrals can be done analytically). Above 
all, they emphasize that DFT should be formulated so that it converges in the limit of a full basis set0 
as obtained for wave function methods (cf. the partial hierarchy of CC approximations: MBPT(2) < 
CCSD < CCSD(T) < full CI). They contrast this precept with conventional DFT approaches such 
as LDA or GGA or hybrid DFT (HDFT, which mixes semi-local GGA and non-local HF exchange 
functionals), which do not consistently predict better answers when applied in that order, even though 
they are claimed to be better approximations. That is, there is no definite hierarchy. 

The Bartlett construction follows most cleanly by requiring that corrections to the density beyond 
the Kohn-Sham density, as calculated according to MBPT, must vanish. As noted in Section 12. 4^ this 
is an alternate but equivalent prescription for orbital-based ab initio DFT. The Hamiltonian is split as 
in the last section into H = Hq + Hi using the Kohn-Sham potential to define the reference state. We 
will use |$Ks) instead of |$o) to make this clear. The perturbative expansion for the energy is [58l [59] 

E = Eo + E^^^ + ^(2) + . . . (31) 
where (in quantum chemistry coupled cluster notation) 

^0 = ($Ks|i^o|$Ks) = , (82) 

i 

E« = (<l>Ks|i^i|$Ks) , (83) 

= (<|.Ks|^1^0^l|$Ks) , (84) 

and so on. Here the operator _Ro is 

^0 = Q , (85) 

Eq — -f^o 

with 

Q = Y,mm + m')m'\ + \^tk)m't\ + --- , (86) 

which are singles, doubles, and higher excitations with respect to the reference KS state (which is a 
determinant). By convention, a, 6, c are unoccupied and are occupied in the KS state. We see 

that this agrees with the Goldstone expansion from before, because Q ensures that only connected 
contributions are included (i.e., it is a complete set of states excluding |$ks))- 

Rather than use the functional derivative to calculate wks, the condition that the density be un- 
changed from the Kohn-Sham density is used. For calculating the exchange-only KS potential \4, this 
implies the density be unchanged to first order in perturbation theory. We expand the wave function 
and the density: 

1^) = |$j,s) + 1^(1)) + 1^(2)) + ... , (87) 
P = PKS + P^'^ + P^'^ + ■ • • , (88) 

where the density is computed from the expectation value of p(x) (note that Bartlett et al. write this 
using 6, which is the same as p): 

^^^^ = ■ ^ ^ 



^^Thc diagonalization of the Hamiltonian including all possible particle-hole excitations in a given basis set is called 
"full CI" . That solution is an upper bound to the exact energy and reproducing it using that same basis set is considered 
to be an unambiguous measure of success in describing electron correlation [77^ . 



26 



So at leading order, we determine fx(x) such that p*^^^(x) = 0. Thus 

($Ks|p(x)|*(^))+c.c. = 0, (90) 
where writing H\'^) = Elip) to this order gives 

{Eo - Ho)\^^''>) = {Hi - i?('))|$Ks) (91) 



and therefore 
Substituting for l^^^^), 



|^«) = ^o(i^i-^^'^)|<fKs). (92) 



p(i) = = (<l>Ks|p(x)i?o^i|<fKs) + c.c. , (93) 
and noting that at this order only single excitations contribute, we obtain 

= X^(^|p(x)|a)({aj|^i|ij} - {a\Vii^\i)) I {ei - Ea) + c.c. 

= -^0*(x)0,(x)((a|^K) + (a|V;K))/(£i-£a) + c.c. . (94) 

i,a 

where we have used ($^|£'o "~ Hq\(^1) = Si — Sa- Here K is the conventional exchange (Fock) energy 
operator VP12, where P12 is the particle-exchange operator. 

Note that we cannot choose {a\K\i) = — (a|\4|i) because must be local. We can think of Eq. (IMI) 
as the solution to a weighted least squares problem to replace the non-local K by the local Vx in the 
space spanned by {(p*, (pa} |5SJ. This is a point-wise identity that can be written [with Xs from Eq. PH]) ] 



j rfx'xs(x,xOK(x') + K(x')] = , (95) 

which agrees with the OEP equation in Section [21 that is, we have rederived the OEP equation with 
exchange only. Here 

/s:(x)0j(x) = y"rfx'p(x,x')l^(x,x')0j(x') (96) 

for a local potential with 

p(x,x') = J]n,0*(x')0.(x). (97) 

i 

Details on solving Eq. 0951) in matrix form are given in Refs. [SHI EH] with pointers to the literature. 

At second order, where the correlation potential first enters, one encounters issues with slow con- 
vergence. Here we simply sketch the problem and proposed solution to give the flavor. Working to 
second-order MBPT, the second order energy is straightforward to write but naturally more complex. 
One insists that p^^^ + p^^^ = and this defines Vc^^ consistent to this order. The problem is that the 
perturbation H[^^ - E^^) = V -{Vh + vP) - E^^^ (note: no K at this order) can have large, diagonal 
contributions 



($Ksip(x)i$n(*ri(^o - H,)-'m)m\H^''^ - e'^'^^'^ m 

{aj\V{l - Pi2)\M) - 6,,{a\K + K\b) - 6at{j\K + . (99) 

Because we can have a = b and i = j, these are diagonal elements of + V^, which make if^^-' a larger 
perturbation than usual [581 El] • The remedy is to resum diagonal and near-diagonal one-particle terms 
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into Hq (this includes a rotation of the KS orbits so that off-diagonal terms in Hq vanish). Most of the 
effect from the resummation comes with the replacement of KS denominator energies Ep {p is either a 
particle or hole) with diagonal Fock matrix elements fpp, where 

fpp = {P\f\p) =^P- {P\K + Kcb) . (100) 

The details are beyond the scope of this review but can be found in Ref. [78] along with results. 

In summary, the key issue is that the freedom in how we split the Hamiltonian can be exploited to 
improve convergence of the perturbation series without sacrificing the form of ah initio DFT. In particu- 
lar, rather than use the standard Kohn-Sham choice (called Gorling-Levy perturbation theory [79| [80]). 
we shift as much physics as possible into Hq. While far from a complete treatment, we hope we have 
conveyed the flavor of the Bartlett approach. As a side note, Bartlett et al. claim that observations 
made for conventional DFT that orbital energies are meaningless (other than the least bound) are 
mostly because functionals and potentials were not accurate. 

3.3 Low- momentum interactions 

The last two sections illustrate how exchange-correlation functions can be constructed for Coulomb 
systems for which many-body perturbation theory is applicable. This is irrelevant for nuclear physics 
unless we have similar control of MBPT. Here we highlight recent results that show that modern RG 
methods can be used to obtain a convergent expansion of nuclear matter properties by evolving to 
low-momentum interactions. There are various methods to achieve these interactions starting from 
phenomenological or EFT potentials (e.g., see Refs. [371 Ell [36]), which for our purposes here are 
equivalent. We simply cite two proof-of-principle results and refer the reader to the literature. 

First we consider the uniform system. Despite a decades-old emphasis on infinite nuclear matter, 
most advances in microscopic nuclear structure theory over the last decade have been through expand- 
ing the reach of few-body calculations. This has unambiguously established the quantitative role of 
three-nucleon forces (3NF) for light nuclei {A ^ 12) [HI [201 [121 [22]. Until recently, few-body fits 
have not sufficiently constrained 3NF contributions at higher density such that nuclear matter calcula- 
tions are predictive. Nuclear matter saturation is very delicate, with the binding energy resulting from 
cancellations of much larger potential and kinetic energy contributions. When a quantitative reproduc- 
tion of empirical saturation properties has been obtained, it was imposed by hand through adjusting 
short-range three-body forces (see, for example, Refs. [83l IM] ). 

Progress for controlled nuclear matter calculations has long been hindered by the difficulty of the 
nuclear many-body problem when conventional nuclear potentials are used. But recent calculations [IQ] 
overcome the hurdles by combining controlled starting Hamiltonians based on chiral effective field theory 
(EFT) [551 [SB] with renormalization group (RG) methods [53, [ST] to soften the short-range repulsion 
and short-range tensor components of the initial chiral interactions [88]. By doing so, the convergence 
of many-body calculations is vastly accelerated [39l [89| [90]. 

Some key results are summarized in Fig. [5l which shows the energy per particle of symmetric matter 
as a function of Fermi momentum /cp, or the density p = 2A;|/(37r^). A grey square representing the 
empirical saturation point is shown in each of the nuclear matter figures. Its boundaries refiect the ranges 
of nuclear matter saturation properties predicted by phenomenological Skyrme energy functionals that 
most accurately reproduce properties of finite nuclei. The calculations of Fig. [5] start from the N^LO 
nucleon-nucleon (NN) potential (EM 500 MeV) of Ref. [85]. This NN potential is RG-evolved to low- 
momentum interactions Viowk with a smooth Wexp = 4 regulator [57]. For each cutoff A, two couplings 
that determine the shorter-range parts of the N^LO 3NF [HH [22] are fit to the binding energy and 
the ^He matter radius [82] using exact Faddeev and Faddeev-Yakubovsky methods as in Ref. [93] . 

The Hartree-Fock results show that nuclear matter is bound even at this simplest level. A calcula- 
tion without approximations should be independent of the cutoffs, so the spread in Fig. [5] sets the scale 
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kp [fm ^] kp [fm ^] kp [fm ^] 

Figure 5: Nuclear matter energy per particle as a function of Fermi momentum kp at the Hartree- 
Fock level (left) and including second-order (middle) and particle-particle-ladder contributions (right), 
based on evolved N^LO NN potentials and 3NF fit to Es^ and r4He- Upper bounds to the theoretical 
uncertainties are estimated by the NN (lines) and 3N (band) cutoff variations [KI] . 

for omitted many-body contributions (more precisely, it sets a lower bound to omitted contributions). 
The second-order results show a dramatic narrowing of this spread, with predicted saturation consis- 
tent with the empirical range. The narrowing happens across the full density range. This is strong 
evidence that these encouraging results are not fortuitous. The particle-particle-ladder sum is little 
changed from second order except at the lowest densities shown. The latter is not surprising because 
at very low density the presence of a two-body bound state necessitates a nonperturbative summation. 
Similar results are obtained using flow equations to evolve Hamiltonians, which is called the Similarity 
Renormahzation Group (SRG) [9ll [951 ESI EZ] in nuclear physics. 

The decrease in cutoff dependence in Fig. |5] with more complete approximations is necessary but 
not sufficient to conclude that the calculations are under control. Indeed, approximations that are 
independent of the cutoff will shift the answer but not widen the error band from cutoff variation. The 
theoretical errors arise from truncations in the initial chiral EFT Hamiltonian, the approximation of the 
3NF, and the many-body approximations. The 3NF approximation is particularly uncertain because 
it involves long-range contributions independent of the cutoff. Many-body corrections to the current 
approximations include higher-order terms in the hole-line expansion and particle-hole corrections. An 
approach such as coupled cluster theory that can perform a high-level resummation including long-range 
correlations will ultimately be necessary for a robust validation. 

These results, while not conclusive, open the door to ab-initio density functional theory (DFT) both 
directly (as in the last sections) but also based on expanding about nuclear matter [511 (next section). 
This is analogous to the application of DFT in quantum chemistry and condensed matter starting with 
the uniform electron gas in local-density approximations and adding constrained derivative corrections. 
Phenomenological energy functionals (such as Skyrme) for nuclei have impressive successes but lack a 
(quantitative) microscopic foundation based on nuclear forces and seem to have reached the limits of 
improvement with the current form of functionals ^71 [HB]- On the other hand, the theoretical errors 
of the calculations in Fig. \5[ while impressively small on the scale of the potential energy per particle, 
are far too large to be quantitatively competitive with existing functionals. However, there is the 
possibility of fine tuning to heavy nuclei, of using EFT/RG to guide next-generation functional forms, 
and of benchmarking with ab-initio methods for low-momentum interactions. Overall, these results are 
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quite promising for a unified description of all nuclei and nuclear matter but much work is left to be 
done. 

Calculations of finite nuclei by another group using soft potentials further support the use of MBPT. 
Roth and collaborators have developed a method using designed unitary transformations to remove 
short-range "hard core" and tensor correlations [99]. The approach is called the "Unitary Correlation 
Operator Method" or UCOM. The result is a soft Hamiltonian that shares the favorable features of 
the RG-based potentials. They have also shown close parallels (and some distinctions) of UCOM to 
the SRG approach |100l llOlj . Using a UCOM NN potential, which is phase equivalent to a given 
initial Hamiltonian (in this case Argonne fig [IB]), they have calculated many nuclei (including heavy 
nuclei such as ^°^Pb) in HF plus second- and third-order MBPT corrections. The energies and radii 
obtained seem to be remarkably converged at second order with good systematics [38]. Preliminary 
results from supplementing the NN potential by a contact 3-body interaction with fitted strength are 
very encouraging [102]. They have also examined corrections from particle- hole states in the RPA for 
closed-shell nuclei; these are found to be relatively small corrections |103j . 

Thus MBPT appears to be a viable candidate for calculating a nuclear exchange-correlation func- 
tional. However, there are three important caveats. Detailed comparisons with a method that does 
high-order resummations will be necessary before one can make robust conclusions. In this regard, 
converged calculations using soft potentials with coupled cluster and no-core-shell-model techniques are 
feasible now in ^^O including 3-body forces and sooioin ^''Ca (NN-only calculations for CC are possible 
now even for an unevolved chiral EFT potential [2S])- Second, tests of MBPT have so far been largely 
limited to spherical nuclei where pairing plays a small role. The final point to emphasize again is that 
corrections are relatively small, but not absolutely small and are large compared to the accuracy sought 
for functionals, so additional contributions must eventually be considered. 



3.4 Density matrix expansion 

While formally the construction of an ab initio DFT energy functional is well-defined based on a many- 
body perturbation expansion (resummed as needed) about a Kohn-Sham reference state, in practice 
the feasibility of using KS potentials from such a functional beyond Hartree-Fock has only recently 
been demonstrated. This is an ongoing area of research in quantum chemistry, with mixed results, and 
progress in nuclear physics will require many further developments. A more immediate route to a DFT 
functional based on microscopic nuclear interactions is to make a quasi-local expansion of the energy 
in terms of various densities, so that functional derivatives needed to define Kohn-Sham potentials 
are immediate. Of course, in doing so one sacrifices the full non-locality present in an orbital-based 
functional. An example of such an expansion is the density matrix expansion (DME) introduced by 
Negele and Vautherin [1041 1105] , which we describe here in some detail. 

The strategy is to follow a path that will be compatible with current nuclear DFT technology but 
testable and systematically improvable. In this regard, the phenomenological nuclear energy density 
functionals of the Skyrme form are a good starting point to build on the MBPT with low-momentum 
interactions. Modern Skyrme functionals have been applied over a very wide range of nuclei, with 
quantitative success in reproducing properties of nuclear ground states and low-lying excitations [5U| 
11061 [TT] . Nevertheless, a significant reduction of the current global and local errors is a major goal [H]. 
One strategy is to improve the functional itself; the form of the basic Skyrme functional in use is very 



^■^These calculations may require an approximation based on a normal-ordering truncation to treat the 3-body force [22j . 
which has not been fully validated for these nuclei. 



30 



restricted, consisting of a sum of local powers of various nuclear densities, e.g., for N = Z nuclei |47] 



where the density p, the kinetic density r, and the spin-orbit density J are expressed as sums over 
single-particle orbitals. (Expressions for the Skyrme functional including isovector and more general 
densities can be found in Ref. [107] .) Fits to measured nuclear data have given to date only limited 
constraints on possible density and isospin dependencies and on the form of the spin-orbit interaction. 
Even qualitative insight into these properties from realistic microscopic calculations could be beneficial 
in improving the effectiveness of the energy density functional. 

A theoretical connection of the Skyrme functional to free-space NN interactions was made long ago 
by Negele and Vautherin using the density matrix expansion (DME) |104l 11051 1108| . but there have 
been few subsequent microscopic developments. The DME originated as an expansion of the Hartree- 
Fock energy constructed using the nucleon-nucleon G matrix [1041 1105] , which was treated in a local 
(i.e., diagonal in coordinate representation) approximation. Recently the DME has been revisited for 
spin-saturated nuclei using non-local low-momentum interactions in momentum representation [96] , for 
which G matrix summations are not needed because of the softening of the interaction (see Section [3^ . 
When applied to a Hartree-Fock energy functional, the DME yields an energy functional in the form of 
a generalized Skyrme functional that is compatible with existing codes, by replacing Skyrme coefficients 
with density-dependent functions. As in the original application, a key feature of the DME is that it is 
not a pure short-distance expansion but includes resummations that treat long-range pion interactions 
correctly in a uniform system. However, we also caution that the Negele- Vautherin DME involves 
prescriptions for the resummations without a corresponding power counting to justify them. 

In essence, the DME maps the orbital-dependent expressions for contributions to the interaction 
energy Ei^^t of the type in Fig. [6](a) into a semi-local form, with explicit dependence on the local densities 
p(R), t(R), V^p(R), and so on. This greatly simplifies the determination of the Kohn-Sham potential 
because the functional derivatives determining the KS potentials can be evaluated directly. The density 
matrix expansion (DME) for spin-saturated nuclei has been formulated for low-momentum interactions 
and applied to a Hartree-Fock energy functional including both NN and NNN potentials in Ref. [96] . 
The output is a set of functions of density that can replace density-independent parameters in standard 
Skyrme Hartree-Fock codes. Furthermore, the upgrade from Skyrme energy functional to DME energy 
functional can be carried out in stages. For example, the spin-orbit part and pairing can be kept in 
Skyrme form with the rest given by the DME. A further upgrade to orbital-based methods would only 
modify the same part of the code, although the increased computational load will be significant. 

Here we outline how the density matrix expansion is carried out for a microscopic DFT at HF order 
using low-momentum (and non-local) two-body potentials. The relevant object we need to expand is 
i?int, which is expressed in terms of the Kohn-Sham orbitals and eigenvalues that comprise the Kohn- 
Sham single-particle propagators as in Section 13.11 For Hartree-Fock contributions only the orbitals 
enter. Higher-order contributions such as the ladder diagrams in the particle-particle channel can also 
be put approximately into the required form by averaging over the state dependence arising from the 
intermediate-state energy denominators (cf. the KLI approximation). Therefore, results presented here 
for the Hartree-Fock contributions to the functional can be generalized to approximately include selected 
higher-order contributions. However, a framework for systematic improvement is yet to be developed. 

Before considering the DME derivation and its application to non-local low-momentum interactions, 
it is useful to first derive in some detail the starting expression for E-^p, the Hartree-Fock contribution 
to i^int. This introduces the basic notation and highlights the differences between most existing DME 





(101) 
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studies, which are formulated with local interactions and in coordinate space throughout, and an ap- 
proach formulated in momentum space and geared towards non-local potentials. For a local potential, 
the distinction between the direct (Hartree) and exchange (Fock) contributions is significant, and is 
reflected in the conventional decomposition of the DFT energy functional for Coulomb systems, which 
separates out the Hartree piece. For a non-local potential, the distinction is blurred because the Hartree 
contribution now involves the density matrix (as opposed to the density) and it is not useful to make 
this separation when the range of the interaction is comparable to the non- locality^ Consequently, 
throughout this section we work instead with an antisymmetrized interaction. 

For a general (i.e., non-local) free-space two-body potential \4[n, -Erf is defined in terms of Kohn- 
Sham states [Eq. ffTTl) ] labeled by i and j, 

A A 

EnF = lJ^{^J\VM^-P^2)\^J) = lJ^{^m^J) . (102) 

ij ij 

The summation is over the occupied states and the antisymmetrized interaction V = Vnn(1 — -P12) has 
been introduced, with the exchange operator P12 equal to the product of operators for spin, isospin, 
and space exchange, P12 = PaPrPr- Note that the dependence of E^p on the Kohn-Sham potential is 
implicit. By making repeated use of the completeness relation in space (r), spin (a), and isospin (r), 

1 = ^ / rfr|rar)(rar| , (103) 

E'hf can be written in terms of the coordinate space Kohn-Sham orbit als as 

Ehf = ^ XI / '^'^W ^"^2 / drs / c/r4 (rio-irir2cr2T-2|V|r30-373r4cr4r4) 

^ ij Wr}-^ J J J 

X 0*(rio-iri)0i(r3O-3r3)0*(r2a2r2)</)j(r4a4r4) . (104) 
From the definition of the Kohn-Sham density matrix, 

A 

P(r3cr3r3, riaiTi) = ^ 0* (ricriri)0i(r3(T3r3) , (105) 

i 

SO Eq. fll04p can be written as 

^^HF = ]^^jdYi---jdri (ricririr2(T2r2|V|r3cr3r3r4a4r4)p(r3a3r3, ricriri)p(r4(T4r4, r2a2r4) 

{o-r} 

= ^TT.Ti^Jdr,--- y"t^r4(rir2|V^®>3r4)p^'nr3,ri)p('Hr4,r2) , (106) 

where a matrix notation is used in the second equation and the traces denote summations over the spin 
and isospin indices for "particle 1" and "particle 2". Hereafter we drop the superscripts on V and p 
that indicate which space they act in as it will be clear from the context. 
Expanding the p matrices on Pauli spin and isospin matrices we have 

P(ri, r2) = ^[po(ri, r2) + Pi(ri, r2)T^ + 5^o(ri, ^2) ■ a + Si{ri, V2) ■ ar^] , (107) 

where we have assumed the absence of charge-mixing in the single-particle states and the components 
are obtained by taking the relevant traces of p. From now on we consider only terms in the energy 
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Figure 6: (a) Schematic diagram for approximations to E^^t that can be expanded using the DME. 
(b) Coordinates appropriate for the DME apphed to the Hartree-Fock potential energy with a non-local 
potential. 

functional arising from products of the scalar-isoscalar (po) density matrices in Eq. (I106p . which are 
the relevant terms for spin-saturated systems with N = Z. Because there will be no confusion, we will 
drop the subscript "0" on the density matrices. 

The first step is to switch to relative/center-of-mass (COM) coordinates (see Fig. [6]). The free-space 
two-nucleon potential is diagonal in the COM coordinate, so the starting point for the DME of the 
two-body Hartree-Fock contribution from a non-local interaction is 

i?HF = ^y dRdrdr' p{R+^-,R+^-)p{R-^-,R-^-)TT,r[{r\V\r')] , (108) 

where V denotes the antisymmetrized interaction and the trace is defined as 

Tr<,,[(r|V|r')] = ^(rairi(T2r2|i>(l - Pi2)Wa,na2T2) . (109) 

{ar} 

The DME derivation of Negele and Vautherin (NV) ^104j focuses on applications to local potentials, 
which satisfy (r|\^|r') = 6{r — r'){r\V\r') . While the original NV work included coordinate-space 
formulas applicable for non-local interactions, for low-momentum potentials it is advantageous to revisit 
and extend the original derivation to a momentum-space formulation. (Note that Kaiser et al. have 
shown how to use medium-insertions in momentum space in their application of the DME to chiral 
perturbation theory at finite density |109l IllOl Ullj .) 

For the momentum space formulation, we first rewrite the density matrices in Eq. (llOSp as 

p(R ± r'/2, R ± r/2) = p(R± ± A/2, R± t A/2) , (110) 

where the vectors appearing on the right-hand side are defined by (see Fig. [6]) 

R± = R±ls, S = i(r' + r), A = l(r'-r). (Ill) 

Introducing the Fourier transform of V in the momentum transfers conjugate to S and A, 

q = k-k', p = k + k', (112) 
(where k', k correspond to relative momenta) gives 

= ^/^^/ ^^(R'q'P)Tw[V(q,p)] , 

^''However, it is useful to separate out the long-distance part of the potential, which is local, and treat its direct 
(Hartree) contribution exactly. 
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where we have defined 



F(R, q, p) = J ciE dA e'"^'^ e'^'^ p(R+ - A/2, R+ + A/2)p(R- + A/2, R" - A 



/2) , (113) 



and 



V(q, p) = 8 j dlldA e"^^-^ e-*P-^ (S - A| V|S + A) . (114) 

The momenta q and p correspond to the momentum transfers for a local interaction in the direct and 
exchange channels. That is, the direct matrix element is a function of q and the exchange is a function 
of p. In contrast, for a non-local interaction the direct and exchange matrix elements depend on both 
q and p. This is the reason why we do not attempt to separate out the Hartree (direct) and Fock 
(exchange) contributions to -Erf, as is commonly done for local interactions. 

The trace of Eq. flll4p can be written in a more convenient form for the application of low momentum 
NN interactions as a sum over partial wave matrix elements (see Ref. [96] for more complete definitions), 

Tr<,,[V(q,p)] = 87r5^ ' {2j + l){2t + I) Pi{^ ■9){klsjt\V\k'lsjt) , (115) 

Isj 

where the primed summation means that it is restricted to values where l + s+t is odd, with k = |(p+q) 
and k' = |(p — q). For simplicity we assume a charge- independent two-nucleon interaction, although 
charge-dependence can easily be included. 

The expression Eq. (11131) for E'hf is written in terms of off-diagonal density matrices constructed from 
the Kohn-Sham orbitals and so is an implicit functional of the density. To circumvent the application 
of the chain rule for the KS potentials, we apply Negele and Vautherin's DME to -Ehf, resulting in an 
expression of the form 

^dme[p, r, J] = j c/R^dme(p(R), r(R), J(R)) , (116) 
with explicit dependence on the local quantities p(R), t(R), and |Vp(R)p that we write for i?HF as 

Ehf = jdn{A[p] + B[p]T + C[p]{Vpf + ■■■). (117) 

(We have suppressed terms that go beyond the present limited discussion; also, when N ^ Z these are 
functions of the isovector densities as well.) The goal is to find the coefficient functions in Eq. (I117p . 
The starting point of the DME is the formal identity |1U4] 

p(R + s/2, R - s/2) = J] r (R + s/2)</.(R - s/2) = [e-(v-v.)/2 ^ ct^^i^Mi^,)]^^^^^^^ , (118) 

a a 

where Vi and V2 act on Ri and Ri, respectively, and the result is evaluated at Ri = Ri = R. We 
assume here that time-reversed orbitals are filled pairwise, so that the linear term of the exponential 
expansion vanishes. After applying a Bessel-function expansion (which is simply the usual plane-wave 
expansion with real arguments), the angle-averaged density matrix takes the form 



p(R + s/2,R-s/2) 



sA;f(R) 



X^(4n + 3)j2n+i(sA;F(R))Qn 



,n=0 



Vi- V2 ' ' 
2A;f(R) 



p(Ri,R2), (119) 



where Q is related to the usual Legendre polynomial by Qi^z^) = P2k+i{iz) / {iz) and an arbitrary 
momentum scale /cf(R) has been introduced. Equation (I119P is independent of k-p if all terms are kept. 
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but any truncation will give results depending on the particular choice for kp. The standard LDA choice 
of Negele and Vautherin, 

fcF(R) = (37rV(R)/2)^/=' , (120) 

is used here. Alternative choices for A;f(R) may better optimize the convergence of truncated expansions 
of Eq. flll9p and lead to a systematic power counting. 

Following Negele and Vautherin, Eq. flllOp is truncated to terms with n ^ 1, which yields the 
fundamental equation of the DME, 

p(R + ^, R - ^) ^ Psl{M^)s) p(R) + s'g{MR)s) [\v'p{R) - r(R) + hp{Ryp{R)] , (121) 
where 

Psl(x) = 3ji(x)/x , gix) = 35Ux)/2x^ , (122) 

and the kinetic energy density is t(R) = |V(/>i(R)p. If a short-range interaction is folded with the 
density matrix, then a truncated Taylor series expansion of Eq. f ll2ip in powers of s would be justified 
and would produce a quasi-local functional. But the local kp in the interior of a nucleus is typically 
greater than the pion mass m^r, so such an expansion would give a poor representation of the physics 
of the long-range pion exchange interaction. 

Instead, the DME is constructed as an expansion about the exact nuclear matter density matrix. 
Thus, Eq. f ll2ip has the important feature that it reduces to the density matrix in the homogeneous 
nuclear matter limit, pnm(R + s/2,R — s/2) = psh{kps) p. As a result, the resummed expansion in 
Eq. fll2ip does not distort the finite range physics, as the long-range one-pion-exchange contribution 
to nuclear matter is exactly reproduced and the finite-range physics is encoded as non-trivial (e.g., 
non- monomial) density dependence in the resulting functional. The small parameters justifying this 
expansion emerge in the functionals as integrals over the inhomogeneities of the density. (See Section 
and Ref. [112] for examples of estimated contributions to a functional for a model problem.) 

In the case of a local interaction, the Fock term is schematically given by Wp ~ J dRds p'^{R + 
s/2, R — s/2)y(s), so a single application of Eq. (I12ip is sufficient to cast E^p into the desired form. For 
a non-local interaction the calculation is more involved as two applications of the DME are required. 
There is arbitrariness in what is kept in higher orders of S^, which is used (following Negele and 
Vautherin) to "reverse engineer" the expansion so that the exact nuclear matter limit is always exactly 
reproduced by the leading term [104] . We emphasize that this is a prescription without established 
power counting or error estimates; as shown in Ref. [9^, different prescriptions can lead to significant 
changes in nuclear observables. The end result for the product of two density matrices is [96| 

p(R + ^, R + l)p{R - ^' - |) ^ pkikFA)p' + ^J:'g{kpJ:) [pV^p psi.{kpA)MkpA) 

- I Vpnjo'(A:FA) +Ji2(A;fA)]) + 2A2^(A;fA)psl(A:fA) QpVV-pr+ ^fc^^^j ^^23) 

In the momentum space expression for -Erf, it remains to evaluate the Fourier transforms defined in 
Eq. flll3p for the expanded density matrices in Eq. fll23p . Identifying the terms in Eq. f lll7p that give 
the DME functionals A[p], B[p], and C[p], we have 

F 

B[p] = -^J^ p'dpTT^r[V{0,p)]l2{p) , (125) 

where Tro-r[V(0, p)] is given by a simple sum of diagonal matrix elements in the different partial waves, 

Tr^.[V(0,p)] = Svr^ '(2j + l)(2t + 1) {llsjt\V\^lsjt) . (126) 

Isj 
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The primed sum is over all channels for which / + s + t is odd. 

In these expressions, the functions Ij{p) and Ij{q) are simple polynomials (and theta functions) in 
the scaled momenta p = p/kp and q = q/kp: 

/Stt 
dxjoipx) = ^(16 - I2p + p') 9(2 - p) , (127) 

/ 357r 
x" dxjoipx) psl{x) g{x) = -^{p'' - W + 40^^ - 2Ap) d{2 - p) , (128) 

h{q) = x'dxjo{qx)g{x) = - — {5f-3)e{l-q), (129) 



hip) = I x^dxjoipx)Ux)psLix) = ^{2-p)e{2-p) , (130) 



hip) = J x^dxUp-x)[foi^)+jl{x)] = ^{A-p')e{2-p). (131) 

Note that the trivial angular dependence of Eqs. f ll27l) -( fT3T|) is a consequence of the angle averaging 
that is implicit with each application of the DME. 

The contributions to Eb^f that have gradients of the local density take the form 

£;hf|,vpP = JdR (Cv2,VV(R) + qvpHVp(R)n . (132) 

We can perform a partial integration on the V^p terms to cast them into the canonical form proportional 
to only |Vpp; that is, 

i?HF|iv,p = / rfR|Vp(R)p [qvpp - ^C7v2j , (133) 

so that 

C[p] = qvpp - ^Cv^p . (134) 

In practice it is efficient and accurate to calculate the derivative in Eq. (11341) numerically rather than 
analytically. The expressions for C|vp|2 and Cyzp are 



1 



2j„ / „2 



1671^4 JO JO 



q dq P dpisiq) hip)Vaviq,p) , (136) 



p 



327r/up Jq 



P^dpl2iP)TTar[ViO,p) 



+ q'dq I p'dphiq) h{p) Vav(g,p) , (137) 







where Va.v{q,p) is the angle- averaged interaction, 

Va.iq,p) = \j rf(cos^^)Tr^,[V(q,p)] , (138) 
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Figure 7: Contribution to the energy per particle in nuclear matter from the isoscalar coefficient func- 
tions A{p), B{p), and C{p) as a function of the density from the DME applied to the Hartree-Fock 
energy calculated using Vfowfc with A = 2.1 fm~^. The result including the NN interaction alone is 
compared to NN plus NNN interactions for two DME expansions (I and II, see Ref. |96j). 

and V(q, p) is given by Eq. f lllSp . Note that care must be taken in the evaluation of dC^2p/dp if the 
vertex V(q, p) is density-dependent or if the local Fermi momentum is not taken to be kp = (37r^p/2)^/^. 

These calculations are extended to (local) three-body forces with corresponding contributions to A, 
B, and C in Ref. ^6] using the same type of expansions. Some representative numerical results given 
there are reproduced in Fig. [TJ which illustrate the relative size of NN and NNN contributions and the 
truncation uncertainty introduced from different DME prescriptions. It is clear that the decreasing hier- 
archy of many-body contributions is maintained in the individual coefficent functions, but cancellations 
magnify the sensitivity to the three-body part and how the truncation is carried out. 

These results are limited and do not yet touch on many of the most interesting aspects of microscopic 
DFT from low-momentum potentials. Topics that should be explored in the future include: 

• Examine the resolution or scale dependence of the energy functional by evolving the input low- 
momentum potential. There will be dependence on the RG cutoff or flow parameter both from 
omitted physics and from intrinsic scale dependence. Calculations at least to second order are 
needed to separate these dependencies. 

• Examine the isovector part of the functional. Contributions from the more interesting long-range 
(pion) parts of the free-space interactions can be isolated, allowing the derivation of analytic 
expressions for the dominant density dependence of the isovector DME coupling functions. 

• Study the dependence of spin-orbit contributions on NN vs. NNN interactions. This includes the 
isospin dependence as well as overall magnitudes. The NN spin-orbit contributions arise from 
short-range interactions, whereas NNN contributions arise from the long-range two-pion exchange 
interaction. Therefore, we expect to find a rather different density dependence for the two types 
of spin-orbit contributions. 

• Explore the contribution of tensor contributions, which have recently been reconsidered phe- 
nomeno logically [1131 1114j . 

• Understand the scaling of contributions from many-body forces. In particular, how does the 
four-body force (which is known at N^LO in chiral EFT with conventional Weinberg counting) 
contribution at Hartree-Fock level impact the energy functional? 

There are both refinements possible within the DME framework and generalizations that test its appli- 
cability and accuracy. Most immediately, the DME at the Hartree-Fock level can be directly extended 



37 



to approximately include second-order (or full particle-particle ladder) contributions by using averaged 
energies for the energy denominators. 

In extending the first calculations from Ref. [96], the standard DME formalism from Negele and 
Vautherin |104] will also need to be modified. This formalism has problems even beyond the truncation 
errors from different DME prescriptions seen in Fig. U\ the most severe being that it provides an 
extremely poor description of the vector part of the density matrix. While the standard DME is 
reasonable at reproducing the scalar density matrices, even here the errors are sufficiently large that 
the disagreement with a full finite-range Hartree-Fock calculations can reach the MeV per particle 
level. Gebremariam and collaborators have traced both of these problems to an inadequate phase 
space averaging (PSA) used in the previous DME approaches |115j . In the derivation of the DME, one 
incorporates average information about the local momentum distribution into the approximation. The 
Negele- Vautherin DME uses the phase space of infinite nuclear matter to perform this averaging for 
the scalar part (and do not even average the vector part). However, the local momentum distribution 
in finite Fermi systems exhibits two striking differences from that of infinite homogeneous matter. 
First, mean-field calculations of nuclei show that the local momentum distribution exhibits a diffuse 
Fermi surface that is especially pronounced in the nuclear surface [116j . Second, the local momentum 
distribution is found to be anisotropic, with the deformation accentuated in the surface region of the 
finite Fermi system [117j . 

Motivated by previous studies of the Wigner distribution function in nuclei, Gebremariam et al. 
calculate the quadrupole deformation of the local momentum distribution using wave functions for that 
nucleus |115] . There are no free parameters. The improvements are substantial for the vector density 
matrices, typically reducing relative errors in integrated quantities by as much as an order of magnitude 
across many different isotope chains |115j . 

Future tests of the DME will include benchmarks against ab initio methods in the overlap region of 
light-to-medium nuclei. Additional information is obtained from putting the nuclei in external fields, 
which can be added directly to the DFT/DME functional. Work is in progress on comparisons to 
both coupled cluster and full configuration interaction calculations as part of the UNEDF project. A 
key feature is that the same Hamiltonian will be used for the microscopic calculation and the DME 
approximation to the DFT. The freedom to adjust (or turn off) external fields as well as to vary other 
parameters in the Hamiltonian permits detailed evaluations of the approximate functionals. In parallel 
there will be refined nuclear matter calculations; power counting arguments from re-examining the 
Brueckner-Bethe-Goldstone approach in light of low-momentum potentials will provide a framework for 
organizing higher-order contributions. These investigations should provide insight into how the energy 
density functional can be fine tuned for greater accuracy in a manner consistent with power counting 
and EFT principles. 

4 DFT as Legendre transform 

In this section, we turn from MBPT in second quantization formalism to the formulation of DFT using 
path integrals for effective actions of composite operators. This can be an intimidating formalism if 
unfamiliar but offers complementary advantages. For example: 

• Effective actions are the natural theoretical framework for Legendre transforms |118l 11191 1120j . 
which is the underlying basis for DFT. These aspects tend to be hidden in the framework of 
Section [31 

• The path integral construction of DFT is transparent, such as the role and usefulness of additional 
densities/sources. (As for implementation, we rapidly discover the same necessity to evaluate 
functional derivatives that was treated in Section [2l) 
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• Path integral effective actions are particularly suited for symmetry breaking, such as encountered 
with pairing. The renormalization issues in the latter case are manifest rather than hidden. 

• Connections to effective field theory (EFT) and power counting can be more accessible. 

• The path integral formulation puts the DFT construction in a broader perspective, which can 
suggest connections and generalizations not apparent otherwise. For example, there are alternative 
effective actions using auxiliary fields or with a two-particle-irreducible nature. The latter may 
be related to more general EDF constructions as proposed in Ref. |16j . 

• The quantization of gauge theories was greatly facilitated by Faddeev-Popov and BRST methods 
using path integrals; the same techniques offer alternative possibilities for implementing collective 
coordinates to restore symmetries broken by the mean- field organization of DFT [12H 1122] . 

• The path integral formulation can suggest different types of nonperturbative approximations, such 
as expansions, that might be needed for extensions beyond low-order nuclear MBPT (e.g., 
to handle long-range correlations). 

Space does not permit a detailed exposition of the path integral formulation of DFT. Furthermore, in 
the short term the construction for nuclear DFT based on MBPT and the DME is probably cleanest 
in the formalism of Section [3l Therefore we focus on presenting the main ideas through analogy and 
schematic versions of the path integrals, and concentrate the EFT discussion on general issues such as 
power counting for functionals. 

As stressed in Section [T], microscopic DFT follows naturally from calculating the response of a many- 
body system to external sources, as in Green's function methods, only with local, static sources that 
couple to densities rather than fundamental fields. (Time-dependent sources can be used for certain 
excited states.) It is profitable to think in terms of a thermodynamic formulation of DFT, which 
uses the effective action formalism [32] applied to composite operators to construct energy density 
functionals |123l 11241 1125] . The basic plan is to consider the zero temperature limit of the partition 
function Z for the (finite) system of interest in the presence of external sources coupled to various 
quantities of interest (such as the fermion density). We derive energy functionals of these quantities by 
Legendre transformations with respect to the sources [29]. These sources probe, in a variational sense, 
configurations near the ground state. 

The work by Lieb [126] on the Hohenberg-Kohn theorem [27| establishes that the real issue for 
DFT is the existence of the Legendre transform F\p\ of the ground state energy as a functional E\v\ 
of the potential. The details involve sophisticated mathematics (e.g., convex-functional analysis) that 
is not readily accessible; we recommend Ref. [29j by Kutzelnigg as a gateway to the mathematically 
rigorous literature behind DFT in terms of Legendre transformations!^^! Fortunately, for our purposes 
the familiarity of physicists with Legendre transforms in the context of thermodynamics is all we need. 
We highly recommend Ref. [26] as an elementary introduction to DFT that carries the analogy to 
thermodynamics through various simple examples. 

Effective actions are a natural framework to implement Legendre transformations, motivate approx- 
imations not obvious in MBPT, and consider generalizations. One limitation of DFT is the exclusive 
role of local potentials (sources) and densities, where locality is in reference to coordinate space. Kutzel- 
nigg points out that this is in contrast to many-body methods that introduce a finite basis in which 
operators are expanded, for which local operators have no privileged place. In this sense, density matrix 
functional theory, as proposed for nuclei in Ref. [T6], seems more natural [29j. By looking at effective 
actions as a broader context, the limitations and problems of local sources are apparent, but also the 
opportunities for generalizations. 

^^There are important formal details [3D], such as that we need E\v\ to be concave to carry out the transform. 
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4.1 Analogy to Legendre transform in thermodynamics 

We first preview DFT as an effective action |125j by recalling ordinary thermodynamics with N particles 
as temperature T — > 0. The thermodynamic potential is related to the grand canonical partition 
function, with the chemical potential /i acting as a source to change N = (N), 

= -kT In Zifi) and N = - ( ^] . (139) 



TV 



Because Q is convex, is a monotonically increasing function of /x and we can invert to find n{N) and 
apply a Legendre transform to obtain 

F{N) = n(/i(A^)) + fi{N)N . (140) 

This is our (free) energy function of the particle number, which is analogous to the DFT energy func- 
tional of the densityO Indeed, if we generalize to a spatially dependent chemical potential J(x), then 

Z(/i) — > Z[J{x)] and fiN = fij^^^P — > j J(x)7/'V(x) • (141) 

Now Legendre transform from lnZ[J(x)] to r[p(x)], where p = and we have DFT with F simply 

proportional to the energy functional! [Note that J(x) — > fext(x) to match our previous notation.] 

The functional F is a type of effective action |125j . An effective action is generically the Legendre 
transform of a generating functional with an external source (or sources). For DFT, we use a source to 
adjust the density (cf. using an external applied magnetic field to adjust the magnetization in a spin 
system). Consider first the simplest case of a single external source J(x) coupled to the density operator 
'p{x) = ip'^{x)il!{x) in the partition function 

Z[J] = e-'^t'^l ~ T:^e-P{H+Jp) _ j v[^^V[^] e-/[^+-^^''^l , (142) 

for which we can construct a (Euclidean) path integral representation with Lagrangian C [32] • (Note: 
because our treatment is schematic, for convenience we neglect normalization factors and take the 
inverse temperature /? and the volume Q equal to unity in the sequel.) The static density p(x) in the 
presence of J(x) is 

p(x) = (p(x)), = , (143) 

which we invert to find J[p] and then Legendre transform from J to p: 

T[p] = -W[J] + Jd^ J(x)p(x) , (144) 



with 

j^^, 6T[p] 6T[p] 



5p(x) (5p(x) 



= . (145) 

Pgs{x) 



For static p(x), F[p] is proportional to the conventional Hohenberg-Kohn energy functional, which by 
Eq. f ll45p is extremized at the ground state density Pgs(x) (and thermodynamic arguments establish 
that it is a minimum |124j )l^ 



^^Because Wext is typically given rather than eliminated, for a closer analogy we would also define ri^(iV) = F{N) — iiN, 
which depends explicitly on both N and fi. This gives the grand potential when minimized with respect to N [26] . 

-•^^A Minkowski-spacc formulation of the effective action with time-dependent sources leads naturally to an RPA-like 
generalization of DFT that can be used to calculate properties of collective excitations. 
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Consider the partition function in the zero-temperature hmit of a Hamiltonian with time-independent 
source J(x) [120j : 

H{J) = H + J Jip^ip . (146) 
// the ground state is isolated (and bounded from below), 



[|0)(0| + C(e-^(^i-^o))] . (147) 



As /? ^ cx), 2[J] yields the ground state of H(J) with energy Eq{J): 

Eo{J) = hm -^logZ[J] = V[J] . (148) 
Substitute and separate out the pieces: 

E,{J) = {H{J))j = {H)j + j J{iPH)j = {H)j + j Jp{J) . (149) 
Rearranging, the expectation value of H in the ground state generated by J[p] i£f| 

{H)j = Eo{J)- j Jp = ^r[p] . (150) 



Now put it all together: 

= 0. (151) 



It[p] = {H)j'-^Eo and J{x] - ^^P] J^o ST[p] 



Pgs(x) 



So for static p(x), r[p] is proportional to the DFT energy functional -Fhk- Furthermore, the true ground 
state (with J = 0) is a variational minimum, so additional sources should be better than just one source 
coupled to the density (as we will consider below) The simple, universal dependence on a non-zero 
external potential v follows directly in this formalism: 

r.M = w,[j]- jjp = w,=^[j + v]- j[{j + v)-v]p = r,=o[p] + J^p- (152) 

Thus allowing for non-zero v^xt is a trivial modification to r[p]. 
There are a various paths to a DFT effective action: 

1. The auxiliary field (Hubbard-Stratonovich) method [V27\ 1128] : Couple ip'^ip to an auxiliary field 
ip, and eliminate all or part of {ip'^ipy. Add a source term Jip and perform a loop expansion about 
the expectation value = A Kohn-Sham version uses the freedom of the expansion to require 
the density be unchanged at each order. 

2. The inversion method |129[ 11241 1130[ [T3T] yields a systematic Kohn-Sham DFT, based on an order- 
by-order expansion. For example, we can apply the EFT power counting for a dilute system. 

3. Derive the functional with an RG approach |132j . This is briefly discussed in Section El 

We will discuss here the inversion method, which connects most closely to the developments in Section [3l 



^^Thc functionals will change with resolution or field redefinitions; only stationary points are observables. This can be 
seen from Eq. (|150p . where r[p] is not the expectation value of H in an eigenstate unless J — J[pgs]- 
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For the Minkowski-space version of this discussion, see Ref. [119] . 
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4.2 Effective actions for composite operators 

A formal constructive framework for Kohn-Sham DFT based on effective actions of composite opera- 
tors can be carried out using the inversion method |123i 11241 \130\ \125\ \133\ \112\ I134[ 1135] . This is an 
organization of the many-body problem that is based on calculating the response of a finite system to 
external, static sources rather than seeking the many-body wave function. It requires a tractable ex- 
pansion (such as an EFT momentum expansion or many-body perturbation theory) that is controllable 
in the presence of inhomogeneous sources, which act as single-particle potentials. As already noted in 
Section 13.31 this is problematic for conventional internucleon interactions, for which the single-particle 
potential needs to be tuned to enhance the convergence of the hole-line expansion [TH [73] , but is ideally 
suited for low-momentum interactions. Given an expansion, one can construct a free-energy functional 
in the presence of the sources and then Legendre transform order-by-order to the desired functional of 
the densities. Because these are complicated, non-local functionals and we require functional derivatives 
with respect to the densities, whose dependencies are usually only implicit, we will ultimately apply 
the methods of Section [2] to derive OEP equations. 

The essential features of the inversion method can be illustrated without the involved path in- 
tegral formalism by considering the Kohn-Luttinger-Ward (KLW) theorem )136l 1137] . which relates 
the perturbative calculation of diagrams using the finite-temperature Matsubara formalism in the zero- 
temperature limit to the calculation of diagrams using zero-temperature perturbation theory. A demon- 
stration of the KLW theorem using an inversion method for the case of an electron gas presented in 
Ref. [31j can be adapted to any hierarchical expansion. Besides perturbation theory in the interaction, 
this includes an EFT expansion relevant to a natural short- distance interaction (which is an expan- 
sion in the Fermi momentum kp times the effective range parameters a^, r^, and so on |138] ) and 
nonperturbative (in diagrams) 1/A^ expansions [139]. 

The basic plan is to carry out order-by-order the conventional thermodynamic Legendre transfor- 
mation already considered: 

F{N) = n{fx) + fxN , (153) 

with fi{N) obtained by inverting N{fi) = —{dQ/dfi)Tv- We expand each of the quantities in Eq. (11531) 
about the non-interacting system: 

Q{fi) = Qoif^) + ^lif^) + ^2{fi) + ■ ■ ■ , (154) 

H = /io + /"i + /U2 H , (155) 

F{N) = Fo{N) + Fi{N) + F2{N) + ■ ■ ■ , (156) 

where the subscript indicates the order of the expansion. (Note that the subscript is just a counting 
parameter that does not have to correspond to a power series in the expansion parameter; e.g., in the 
Coulomb case the expansion parameter is but Q2 has both an term and the correlation energy 
of order e^lne.) The non-interacting system refers to a system of zeroth order in the EFT expansion 
parameter. This means the zeroth-order system has no internal interactions among the particles, but 
it can include external sources (we exploit this freedom below). 
The number of particles is 

Note that we could simply use the unexpanded first equality in Eq. f ll57p together with the series in 
Eq. (1154p . because they define a parametric relation between and Q in terms of /i [31]. Instead, 
to mimic DFT we perform the inversion in Eq. fll57p order-by-order, treating A^ as order zero in the 
expansion. (That is, we ensure there are no corrections to A^ at higher order.) This means that the 
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zeroth order equation, 



N 



:i58) 



is the only equation to which N contributes (by construction). Thus the "exact" is obtained from 
the non-interacting thermodynamic potential This might not sound impressive, but the analogous 
situation holds when we generalize /i to be position dependent or coupled to the pair density in a finite 
system. In these cases, it is the exact, spatially dependent fermion or pair density (with appropriate 
renormalization conditions, see below) that is obtained from a non-interacting system with a single- 
particle potential that is the generalization of no- This is precisely the description of the Kohn-Sham 
system (see Refs. |133] and |112j for details on carrying out the DFT case without pairing). 

Equation fll58p determines A^(/^o) at any temperature, from which we can find fio{N) for any system 
for which we can identify Qq (the inversion is unique because fiQ is a monotonic function of N [31]). 
If we have a uniform system with no external sources, /iq is the chemical potential of a noninteracting 
Fermi gas at temperature T with density N/V. In particular, at T = with no external potential and 
spin-isospin degeneracy u, 

^o(iV) = {6n^N/uVf/^ = kl/2M = e° . (159) 

The first-order equation extracted from Eq. (I157P has two terms, which lets us solve for /xi in terms of 
known (from diagrams) functions of fiQ-. 
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[dQi/dfi] 



(160) 



At second order, we can isolate and solve for /i2, eliminating fii using fll6Up . This pattern continues to 
all orders: fii is determined by functions of /xq only. 
Now we apply the inversion to F = Q + fiN: 



F{N) = no{fio)+fXoN + /ii 



+ (^i(/io) +/iiAr 



+ yU2 



dflo 

9/i 



1 2 



(161) 



F2 



The fiiN term always cancels with fii[dQo/dfi]^=^Q in Fi for z > 1 because of Eq. 



leaving 



F{N)=Fo{N)+ni{ixo) + nM 

Fi 



1 [dQi/dfi] 



2 [d^no/djj^] 



(162) 



F2 



where we have also used Eq. fll60p to simplify F2. The expansion for F can be extended systematically, 
but this is all we need here. (Higher orders can be found by following the prescription in Refs. 1140(1141] .) 

This construction is rather general. The Kohn-Luttinger-Ward theorem explores a particular case, 
the T — > limit, in which the second term in F2 in Eq. fll62p cancels precisely against the anomalous 
diagram in ^2, as illustrated in Fig. [81 This cancellation of derivative terms and anomalous diagrams 
occurs to all orders in the expansion. An analogous cancellation occurs in the Kohn-Sham DFT inver- 
sion [133j . The end result is an expression for the free-energy F{N) in terms of the diagrams used for 
f2j(yLi), only evaluated with fi = fiQ and excluding the anomalous diagrams (both of which simplify the 
evaluation of F{N)\). This is precisely the formalism used in Ref. |138j for a uniform low-density Fermi 
gas at zero temperature, where fio appeared as the Fermi energy of Eq. (11590 . 
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Figure 8: Cancellation of the anomalous diagram at NLO. The double lines represents the inverse of 
[92^]o/V]/.=Mo in Eq. 

This procedure can be generalized rather directly [1231 1124] to carry out the inversion from p[J] to 
J[p] needed for Eqs. f ll42p - fll45p . The idea again is to expand the relevant quantities in a hierarchy, 
here labeled by a counting parameter A, 

W[J,X\ = Wo[J] + XWi[J] + X^W2[J] + ■ ■ ■ , (163) 
J[p, A] = Jo[p] + \Ji[p] + AV2[p] + ■ ■ ■ , (164) 
T[p, A] = ro[p] + \T,[p] + X^T^lp] + ■ • • , (165) 

treating p as order unity (which is the same as requiring that there are no corrections to the zero-order 
density), and match order by order in A to determine the Jj's and Fj's. Zeroth order is a noninteracting 
system with potential Jo{x): 

To[p] = -Wo[Jo] + J rfx Jo(x)p(x) (166) 

and 

P(x) = -JYTZV ■ ^^^^^ 

Because p appears only at zeroth order, it is always specified from the non-interacting system according 
to Eq. (11671) : there are no corrections at higher order. This is the Kohn-Sham system with the same 
density as the fully interacting system. 

What we have done is use freedom to split J into Jo and J — Jq, which is the analog of introducing 
a single-particle potential U and splitting the Hamiltonian according to H = {T + U) + {V — U), as 
discussed in Section |3l Typically U is chosen to accelerate (or even allow) convergence of a many-body 
expansion (e.g., the Bethe-Brueckner-Goldstone theory [711 [72], [73]). ^oi DFT, we need to choose it to 
ensure that the density is unchanged, order by order. Thus, we need the flexibility in the many-body 
expansion to choose U without seriously degrading the convergence; such freedom is characteristic of 
low-momentum interactions. (Note: If there is a non-zero external potential, it is simply included with 
Jo-) 

The path integral defined by Wq (or Zq actually) is a gaussian, which is equal to a functional 
determinant that is evaluated by diagonalizing PFo[<^o]- To do so means to introduce Kohn-Sham orbitals 
(pi and eigenvalues e,, 

[-VV2m - Jo(x)]0i = Ei^Pi (168) 

so that 

A 

p(x) = 5^|0,(x)r (169) 

i=l 

Then Zq is the product and Wq is the sum of the e^'s. Thus in the path integral formalism, the Kohn- 
Sham system arises naturally as the zeroth-order approximation to the problem. The organization is 
based on a saddlepoint evaluation about the system defined by Jo (which still must be specified) and 
subsequent corrections. 

The orbitals and eigenvalues are used to construct the Kohn-Sham Green's functions, which are 
used as the propagator lines in calculations the diagrams generated by VFj[Jo]. Finally, we find Jq for 
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the ground state by truncating the chain at Tj^^^, 

j^^W,^V,^J,^W2^T,^...^ ^ r,_^ (170) 
and completing the self-consistency loop that enforces J(x) = 0: 

Jo(x) = - )^ J.(x) = L ^ = ■ (171) 

Note that the sum of the Fj's is directly proportional to the desired energy functional. When transform- 
ing from Wi to Fj, there are additional diagrams that take into account the adjustment of the source 
to maintain the same density and also so-called anomalous diagrams (these are two-particle reducible). 
This is the most complicated part but corresponds directly to the extra terms in the KLW inversion 
[see Eqs. (116ip . (11621) and Fig. [8]. A general discussion and Feynman rules for these diagrams are given 
in Refs. [IMl [M [TSU [T33J. 

We emphasize that even though solving for Kohn-Sham orbitals makes the approach look like a 
mean-field Hartree calculation, the approximation to the energy and density is only in the truncation 
of Eq. fll7ip . It is a mean-field formalism in the sense of a conventional loop expansion, which is 
nonperturbative only in the background field while including further correlations perturbatively order- 
by-order in loops. The special feature of DFT is that the saddlepoint evaluation applies the condition 
that there are no corrections to the density. Once again, this is not ordinarily an appropriate expansion 
for internucleon interactions; it is the special features of low-momentum interactions that make them 
suitable. 

To generalize the energy functional to accommodate additional densities such as r and J (which 
appear in the density matrix expansion for nuclei), we simply introduce an additional source coupled to 
each density. Thus, to generate a DFT functional of the kinetic-energy density as well as the density, 
add ^^(x) V'?/'^ • V?/^ to the Lagrangian and Legendre transform to an effective action of p and r |134] : 

F[p,r] = W[.lr,] - [ dx J(x)p(x) - / rfxr/(x)r(x) . (172) 



The inversion method results in two Kohn-Sham potentials. 



Jn X 



and rjo{x.) 



5Fint[p, r] 



<5p(x) 

where Fjnt = F — Fq. The Kohn-Sham equation is now [134] 



5r(x) 



5 

P 



(173) 



["^M^^ " = "^'^^ ' ^^^^^ 

with an effective mass l/2M*(x) = 1/2M — ?7o(x), just like in Skyrme HF (see also Ref. [142] for an early 
application to Coulomb systems). Generalizing to the spin-orbit or other densities (including pairing 
1 143] , see Section 15. ip proceeds analogously. We note that the variational nature of the effective action 
implies that adding sources will always improve the effectiveness of the energy functional. 

The Feynman diagrams for Wi will in general include multiple vertex points over which to integrate. 
Further, the dependence on the densities will not be explicit except when we have Hartree terms with 
a local potential (that is, a potential diagonal in coordinate representation). One way to proceed is to 
calculate the Kohn-Sham potentials using a functional chain rule, e.g.. 
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Figure 9: Schematic representation of Eq. fll75p for a local potential, where the double- line symbol 
denotes the {6p/6Jo)~^ term. 



and steepest descent |124] . This is illustrated schematically for a local interaction in Fig. [H We see 
that the Kohn-Sham potential is always just a function of R but that the functional is generally very 
non-local. If zero-range interactions are used, these diagrams collapse into an expression for Jo(R) 
that has no internal vertices, but this is no longer true for diagrams with more than one interaction. 
The orbital-based methods from Section [2] take the chain rule in Eq. fll75p one step further, adding a 
functional derivative of the sources with respect to the 0j's (and e^'s) [H EH [5ll 1144] . which can be 
applied directly to the functionals. 



4.3 EFT and power counting for functionals 



The expansion suggested by low-momentum interactions is perturbation theory in powers of the softened 
interaction V. But is is also instructive to consider the simplest EFT example, a dilute system of 
fermions in a harmonic trap, interacting via natural-sized contact interactions [133illl2]ll34j . (Natural- 
sized means that the scattering length is not fine-tuned to a large value compared to the effective range.) 
We can construct the effective action as a path integral by finding W[J] order- by-order in an EFT 
expansion. For a dilute short-range system, the EFT Lagrangian for a nonrelativistic spin-1/2 fermion 
field with spin- independent interactions in Euclidean form is (t^ is the Euclidean time) 
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dtE 
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(176) 



where V= V — V is the Galilean invariant derivative and h.c. denotes the Hermitian conjugateo The 
terms proportional to C2 and contribute to s-wave and p-wave scattering, respectively, while the 
dots represent terms with more derivatives and/or more fields, as well as renormalization counterterms. 



The Euclidean generating functional with chemical potential /i and external sources ri(x) and rj' 



x 



[821 [128]: 



■ fi tpt {x)tpc, {x)+ril, {x)4)a {x)r)a [x)] 



(177) 



where J df"x includes a dtE integration that runs from —{3/2 to (3/2 (to facilitate the (3 ^ 00 limit) and 
anti-periodic boundary conditions are imposed. A conventional perturbative expansion is realized by 
removing the interaction terms from the path integral in fll76p in favor of functional derivatives with 
respect to r] and and performing the remaining Gaussian integration over and tp'^ [32l 1128] : 



(178) 
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We use units with h = I. 
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where the spin indices are imphcit and we have introduced the noninteracting partition function Zq. 
Exphcit expressions for the Green's function in coordinate and momentum space can be found in 
Ref. [32] • The hnked-cluster theorem [32] shows that the difference of the interacting and noninteracting 
thermodynamic potentials Q and Qq is given by the sum of connected diagrams from the expansion of 
Z, with the external sources rj'^ and t] set to zero at the end: 

T, V) - ^]o(/i, T, V) = ^iW[0, 0; /i] - Wo[0, 0; ^]) . (179) 

In evaluating the Feynman diagrams for W[J] and new diagrams for r[p] order by order in the 
expansion (e.g., EFT power counting), the source Jo(x) is now the background field. This means that 
propagators (lines) in the background field Jo(x) are 



G°s(x,x» = 5;]0,(x)0*(x') 

i 

where dAx) satisfies: 



{£i - Sf) _^ 0{eF - £i) 



uj — Si + if] uj — Ei — irj 



;i80) 



+ Wext(x) - Jo(x)] 0i(x) = ei0i(x) . (181) 



L 2M 

For example, if we apply this prescription to the short-range LO contribution (i.e., Hartree-Fock), we 
obtain 

WAM = -u{u-l)cjd'^ ^ |^G^s(x,x;c.)G°,s(x,x;c.O 

^^""'^^C, /'d3xK(x)f , (182) 



2 V 

where v is the spin-isospin degeneracy and 



p,„(x)^z.5^|0.(x)r (183) 

i 

Expressions for the other Wj's proceed directly from conventional Feynman rules using the new propa- 
gator. 

Given W\J\ as an EFT expansion, we perform the Legendre transformation, 

r[p] = W\J\ - Jjp, (184) 

by using the EFT power counting with the inversion method as described above, which gives us the 
means to invert to find J[p] and to make an order- by-order inversion from W[J] to r[p] |129l 11241 1130j . 
It proceeds by decomposing J(x) = Jo(x) + Jlo(x) + Jnlo(x) + . . . as described ear her and fixing Jq 
using 

P^^^ = TTTY ^o(x)U^^ = %ffM . (185) 

Evaluating the functional derivatives is immediate if F is approximated so that the dependence on the 
density is explicit, as with the LDA or DME (see below). Otherwise we apply the OFF chain rule. 

Consider the T = local density approximation (LDA) for a dilute fermion gas with natural short- 
ranged interactions (meaning the scattering length is comparable in magnitude to the effective range 
tq and the p-wave scattering length Op |138] ). The energy density in the uniform system evaluates to: 



E k 



V ^2M 



^ + - ^)§^ikFao) + iiy- 1)^(11 - 2 ln2)(A;Fao)' (186) 



+ (z/ - 1) (0.076 + 0.057(z/ - 3)) {kpaof + (z^ - l)-^(fcFro)(fcFao)' + (z^ + l)-^(/cFap)^ + ■ 

lOvr oTT 
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where kpa^, kpro, and kpUp are expansion parameters with kp = (Gvr^p/i/)^/^. Using this relation to 
replace kp everywhere by p(x), we directly obtain the LDA expression for T[p\, 



r[p] 



M 



+ ih 4[l,(^)f'' + 4 4 r„lp(x)]»'^' + <i»1p(x)]»'^' + 



8/3 



18/3 



;i87) 



The Kohn-Sham Jq according to the EFT expansion follows immediately in the LDA from fllSSp : 



^o(x) 



1) A-nar 
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M 



P(x) 



Cl 



2M 



[p(x)]^/=^ - C2 al\p{^)fr^ - C3 ai ro[p(x)]^/^ - c, a^[p(x)]^/^ + 



15/3 



5/3 



15/3 



where the {(ij}'s and {cj}'s are given in Ref. [145] . Sample results at different EFT orders for a dilute 
Fermi gas in a harmonic oscillator trap are shown in Fig. [TOl Note the systematic convergence with 
successive orders in the LDA average (kptts) (and = Oq) for both the energy and density. 

An important consequence of the systematic EFT approach is that we can also estimate individual 
terms in energy functionals. If we scale contributions to the energy per particle according to the average 
density or (kp), we can make such estimates |133[ I112j . This is shown in Fig. [11] for both the dilute 
trapped fermions, which is under complete control, and for phenomenological energy functionals for 
nuclei, to which a postulated QCD power counting is applied [146[ I147 J. In both cases the estimate^ 
agree well with the actual numbers (sometimes overestimating the contribution because of accidental 
cancellations), which means that truncation errors are understood. Understanding how such power 
counting can emerge from the MBPT-based DFT for nuclei is an important problem for future study. 



4.4 Additional comments 

The perspective of Legendre transforms clarifies many DFT issues |29] . For example: 

• The Legendre transform formulation adds mathematical rigor to DFT claims. For example, one 
finds that a bijective mapping between v and p does not exist in general. However, it is not 
necessary for there to be such a mapping to have a Legendre transform; only the uniqueness of p 
given V is needed. 

^^The symbols with error bars are natural estimates with a factor of two spread in the order unity coefficient. 
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Figure 11: Estimates for energy functionals for a dilute fermions in a harmonic trap (left) and for 
phenomenological energy functionals for nuclei (right). 



• For non-degenerate states, it is true that the density is the functional derivative of the energy 
with respect to the external potential. But this is not actually needed and so is not a problem for 
degenerate states; only the concavity of E[v] is necessary. 

• DFT is only variational with respect to p if we have the exact functional (i.e., never). Once there 
are approximations, it is no longer variational. This is not a cause for concern, as other successful 
many-body approaches such as coupled cluster are also not variational. 

Additional useful comments on DFT and Legrendre transforms can be found in Ref. [29] . 



5 Topics for nuclear DFT 

In this section, we consider various topics specifically relevant to implementations of DFT for nuclei. 
5.1 Pairing 

Pairing is an important feature of finite nuclei that is absent from isolated atoms and molecules. Its 
inclusion in nuclear DFT based on microscopic interactions is a topic of much current activity. While 
there is some work on superconductivity in a Coulomb DFT framework, it is based on using non-local 
source terms to avoid divergences. The use of local densities for pairing is generally preferred for finite 
nuclei from a computational point of view. However, we note that recent results from hybrid EDF 
calculations using low-momentum potentials at lowest order, which suggest that an accurate ab initio 
DFT treatment of pairing is feasible [IS] , use non-local pairing fields. To use local pairing densities, we 
need a consistent generalization of the Skyrme-Hartree-Fock-Bogoliubov approach |lU6j . In fact, one 
can use functionals with local pairing fields as in phenomenological Skyrme functionals (that are based 
on zero-range effective interactions) by properly renormalizing [521 1148] . 

Pairing is an example of spontaneous symmetry breaking, which is naturally accommodated in an 
effective action framework [118^ 1119] . For example, consider testing for zero-field magnetization M in a 
spin system by introducing an external field H to break the rotational symmetry. Legendre transform 
the Helmholtz free energy F{H): 

invert M = -dF{H)/dH =^ T[M] = F[H{M)] + MH{M) . (189) 

Because H = dT/dM — > 0, we look for the stationary points of F to identify possible ground states, 
including whether the symmetry broken state is lowest. For pairing, the broken symmetry is a U{1) 



49 



[phase] symmetry for fermion number. We apply an external source that breaks the number symmetry, 
forcing ipip to have an expectation value. Then we turn that source off and see whether the expectation 
value (condensate) persists. (Note: we will only have actual spontaneous symmetry breaking in an 
infinite system.) 

The textbook effective action treatment of pairing in condensed matter is to introduce a contact 
interaction [12711128] : g ip^i/j^ipip, and perform a Hubbard-Stratonovich transformation with an auxiliary 
pairing field A{x) coupled to ip^ip^ which eliminates the contact interaction. Then one constructs the 
one-particle-irreducible effective action r[A] with A = (A), and looks for values for which 6r/6A = 0. 
To leading order in the loop expansion (which is a mean field approximation), this yields the BCS 
weak-coupling gap equation with gap A. 

The natural alternative here is to use the inversion method for effective actions again [ 129[ll24]ll30j . 
Thus we introduce another external current j(x), which is coupled to the fermion pair density ipip in 
order to explicitly break the phase symmetry. This is a natural generalization of normal-state Kohn- 
Sham DFT [SH 1149] . The generating functional has sources J,j coupled to the corresponding 
densities |143j : 

Z[JJ] = e"^['^'-''l = / D{ip^ip) e'-^'^*'' [C+J(x)i^iip^+j{x){4,\ipl+4,i4,T;)] ^ (190) 



(191) 



Densities are found by functional derivatives with respect to J and j: 

6W[J,j] 



p{x) = {tlj\x)^{x))jj 



6J{x) 

and (note that k is called in Ref. |143] ) 

5W[J,j] 



K 



[x) = {^\{x)^Ij\{x) + ^i{x)^Ij^{x)) j^j 



Sj{x) 



(192) 



(The source j would in general be complex, but it is sufficient for our purposes to take it to be real.) 
The effective action r[p, k] follows as in Section H] by functional Legendre transformation: 

T[p,k] = W[J,j] - J (fxJ{x)p{x) - j d^xj{x)K{x) , (193) 

and is proportional to the (free) energy functional E[p,k]] at finite temperature, the proportionality 
constant is j3. The sources are given by functional derivatives wrt p and k: 

But the sources are zero in the ground state, so we determine the ground-state p(x) and k(x) by 
stationarity: 



6E[p,K] 



6E[p,K] 



. (195) 



5p(x) 

This is Hohenberg-Kohn DFT extended to pairing! 

We need a method to carry out the Legendre transforms to get Kohn-Sham DFT; an obvious choice 

is to apply the Kohn-Sham inversion method again, with order-by-order matching in the counting 
parameter A. Once again, 

diagrams =^ W[J,jA]=Wo[J^j] + \Wi[J,j] + \^W2[J,j] + --- (196) 

assume =^ J\p, k,X\ = Jq[p, k] + \Ji\p, k] + J2[p, k] + ■ ■ ■ (197) 

assume =^ k. A] = jo[p, + Aji[p, k] + A^JsIp, H (198) 

derive ^ F [p, A] = Fq [p, + AFi [p, k] + A^Fs [p, + ■ • • (199) 
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Start with the exact expressions for F and p 



and 



r|p,K] = »'[j,j]- jjp- jjn 



(200) 



(201) 



5 Ax) ' ^ ^ 6jix) ' 

and plug in the expansions, with p, n treated as order unity. Zeroth order is the Kohn-Sham system 
with potentials Jo(x) and io(x), 



To[p,k] = Wo[Jo,jo] - jjop- jjofi 



so the exact densities p(x) and k(x) are by construction 



I . _ 5Wq[JoJg] I ^ _ W[^o, jo] 



6Jo{x) ' Sjo{x) 
Now introduce single-particle orbitals to diagonalize Fq, which means solving 

Ui(x) 



/io(x)-po jo(x) 

jo(x) -/io(x)+/io 



ni(x) 
?^i(x) 



where 



/in fx) = 



2M 



+ Wext(x) - Jo(x) . 



(202) 



(203) 



(204) 



(205) 



As expected, this is just like the Skyrme Hartree-Fock Bogoliubov (HFB) approach [17] . For ab initio 
DFT based on wave-function methods, Bogoliubov transformations lead to the same generalized Kohn- 
Sham equations. 

The diagrammatic expansion of the Wi^s is the same as without pairing, except now lines in diagrams 
are KS Nambu-Gor'kov Green's functions [150] . 



G 



/ {Tij^{x)^Ux'))o {T^^{x)^i{x')), 
[{T^l{x)^\{x'))o {T^l{x)^i{x'))o 



(206) 



where the time-ordering is with respect to Euclidean time. In frequency space, the Kohn-Sham Green's 
functions are 



3 

F°(x,x» = -5^ 



m,-(x)m*(x') ^ i;,-(x')i;*(x) 



iuo — Ej 



iuj — 



iuo + Ej 



iu + E^ 



(207) 



(208) 



The Kohn-Sham self-consistency procedure involves the same iterations as in phenomenological Skyrme 
HF (or relativistic mean-field) when pairing is included. In terms of the orbitals, the fermion density is 



p(x) = 2 5^ |t;,(x) 



(209) 



and the pair density is (warning: this is unrenormalized!) 



(210) 
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Figure 12: (a) Feynman diagram at second order in a perturbative expansion of j] in j. (b) Con- 
verence of the integral for the pair density in the uniform system for two subtractions as a function of 
an energy cutoff in the integraL 



The chemical potential /xq is fixed by J p(x) = A. Diagrams for r[p, k] oc Eq[p, k] + -Eint[p, k\ yield the 
Kohn-sham potentials 



Jo(x) 



P=Ps 



5p(x) 



and jo(x) 



P=pgs 



(211) 



So it appears that everything goes through smoothly as in the ab initio DFT without pairing. 

Potential problems arise, however, because the use of local composite operators leads to new ultra- 
violet divergences even at the mean-field (Hartree-Fock) level when pairing is included. Divergences 
at this order are not encountered when coupling an external source to the fermion density ip'^ip, but 
appear now because the composite operators ipip and ip'^ip'^ need additional renormalization [151] . The 
divergences at leading order are symptomatic of generic problems identified long ago by Banks and Raby 
|152] that arise with effective potentials of local composite operators. (It was to avoid these divergences 
in DFT that non-local sources were used in Refs. [153[ I154j .) These problems inhibited for many years 
the use of effective actions of composite operators but more recently Verschelde et al. |155i I156i 1157] 
and Miransky et al. [158[ 11591 1160] have revived their use for relativistic field theories. In Ref. [143] . 
it was shown how to identify and renormalize up to NLO the EFT from Section 14.31 when pairing was 
included. 

The source of the divergences is found immediately when we try to carry out the DFT pairing 
calculation even for a uniform dilute Fermi system. The generating functional with constant sources p 
and j for the EFT from Section 14.31 is: 



-W[p,j] 



(212) 



(cf. adding an integration over an auxiliary field JD{A*,A) e 



'\Co\ 



then shifting variables to elim- 



inate '?/'|'?/'|'?/'|'?/'| for A*ipiipi). There are new divergences because of j, e.g., expand W to C(j^) [see 
Fig. [T2r a)]. which has the same linear divergence found in 2-to-2 scattering. ( Equivalent ly, in coordinate 
space there is a l/|ri — divergence in 'y*(ri)Mj(r2) as it becomes local with |ri — 0.) To 

renormalize, we add the counterterm |C|jP to C (see |120j ). which is additive to W (cf. |Ap), so there 
is no effect on scattering. 
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We can follow a less formal and more numerically suitable procedure to renormalize the pair (anoma- 
lous) density, 

= 5Z K(x)^»(x) + Mi(x)t;*(x)] — > oo , (213) 

i 

which diverges for contact interactions in a finite system. This is to use the renormalized expression for 
K in the uniform system, 




= / JO I --o\^ finite, (214) 

which is cut off at momentum kc, and apply this in a local density approximation (i.e., Thomas- Fermi) : 
«:(x)=2^«,(x)t;,(x)-jo(x)^^ with = M_Z + j(x) - ;,o • (215) 

i 

This procedure was worked out by Bulgac and collaborators in Refs. [52 | [5T | ri49i 1161] . Convergence is 
very slow as the energy cutoff Ec is increased, so Bulgac and Yu devised a different subtraction. 



(fk . 1 V 



A comparison of convergence in the uniform system for the two subtraction schemes (12141) and (12161) [see 
Fig. [12] shows dramatic improvement for the Bulgac/Yu subtraction. Bulgac et al. have demonstrated 
that this works in finite systems [149j and it has been adopted for some Skyrme HFB implementations. 

5.2 Broken symmetries 

Ordinary nuclei are self-bound, which presents conceptual issues about whether Kohn-Sham DFT is well 
defined and practical problems on how to deal with the consequences of symmetry breaking by the KS 
potentials, which will not have all of the symmetries of the Hamiltonian [16]@ These broken symmetries 
include the U{1) phase symmetry for fermion number and translational and rotational invariance. While 
restoring broken symmetries is a topic well-explored for mean- field approximations [l6l HT] , it has only 
recently been considered in a context relevant to ah initio DFT. Because no proven practical approaches 
are yet available, we simply point out the issues and current references. 

The textbook discussions of how to restore mean-field broken symmetries tend to follow one of these 
two related lines of discussion: 

• States related by a unitary transformation U{a) corresponding to a broken symmetry are degen- 
erate: 

= [/(«)!</.) (217) 

with 10) a "deformed" state, implies 

{(t>a\H^\ct>a) = {(t>\HN\cl>) . (218) 

The degeneracy can be removed by diagonalizing in the subspace spanned by the degenerate 
states. The group parameter a for continuous groups can be considered a collective coordinate, 
which specifies the orientation in gauge space of the deformed state A general strategy is to 
transform from 3^4 particle coordinates into collective and internal coordinates [47j. 



^^A state is one of broken symmetry if it does not have the quantum numbers of the eigenstates of the Hamiltonian 
(parity, particle number, angular momentum, linear momentum, isospin, and so on). Note that this does not mean it 
needs to be invariant, just that it can be labeled by definite quantum numbers |46j . 
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• In finite systems, broken symmetries arise only as a result of approximations. This usually hap- 
pens with variational calculations over trial wavefunctions that are too restricted; a mean-field 
approximation is an example. The symmetry can be restored by using a linear superposition of 
degenerate states: 



which when minimized with respect to the /(a)'s projects states of good symmetry |16]. (Because 
minimizing with respect to |0) and with respect to /(a) do not commute, there are two types of 
projection. It is most accurate to project first and then find the best deformed state corresponding 
to a given quantum number.) For example, particle number projection for EDF's is described in 
Refs. [T62l[T63l . 

When implemented, these approaches are sometimes considered to be beyond EDF, where there are only 
densities and not a wavefunction. From a different perspective, the restoration of broken symmetries 
of GCM-type configuration mixings should be considered as a "multi-reference" extension of the usual 
"single- reference" EDF implementation (see Refs. [I2l HHl Il9] ) . What about ab initio DFT as we have 
considered it? 

For nuclear DFT, the conceptual question was highlighted by J. Engel |164j . who pointed out that 
the ground state of a self-bound system, with a plane wave describing the center of mass, has a density 
distributed uniformly over space. Clearly this is not the density one wants to find from DFT, so there 
is a question of how to proceed. There are two separate considerations: i) Does Kohn-Sham DFT exist 
in a useful form for self-bound systems? ii) If so, how does one formulate and implement it? Engel and 
other authors have addressed this issue recently |164[ 11651 1166[ 11671 1168[ I169j , with a consensus that 
HK existence proofs for DFT are still well founded, but for internal densities (meaning independent 
of the center-of-mass motion when considering broken translational symmetry) o These authors also 
propose various schemes to carry out Kohn-Sham DFT that should be testable in the near future. 

We have considered ab initio DFT from two viewpoints: MBPT using wave-function methods and 
effective actions with path integrals. How are these symmetry issues dealt with in these approaches? 
Wave function methods have several related strategies for dealing with the "center of mass" (COM) 
problem: 

1. Isolate the "internal" degrees of freedom, typically by introducing Jacobi coordinates. Then the 
observables are by construction independent of the COM. This gets increasingly cumbersome with 
greater numbers of particles. 

2. Work in a harmonic-oscillator Slater determinant basis, for which the COM decouples, and intro- 
duce a potential for the COM that allows its contribution to be subtracted. 

3. Work with the internal Hamiltonian (i.e., subtract the COM kinetic energy Tqm) so that the 
COM part factors and does not contribute to observables to good approximation (see in particular 
Ref. |170j for coupled cluster calculations). 

Versions of the first two possibilities are in fact among the ideas considered for DFT in Refs. |164[ I165[ 



For the effective action approach, the issue of broken symmetry is familiar from the study of soli- 
tons |17H [32] . where it arises as the problem of dealing with zero modes when calculating quantum 
fluctuations. The methods found in the literature are similar to the textbook approaches cited above. 
One compelling approach for ab initio DFT in the effective action formalism uses the Fadeev-Popov 

^•^In other contexts, such densities are cahed "intrinsic", but this has a different meaning in the context of symmetry 
breaking, so "internal" is typicaUy used instead. 




(219) 



[1661 [T671 [T681 [T69 ; . 
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construction (or BRST invariance |121] ) to introduce collective coordinates through ghost degrees of 
freedom. This is worked out to some degree by Calzetta and Hu for broken particle number symmetry 
in Ref. |122j . but a concrete implementation appropriate for nuclei remains to be formulated. 

Thus, while dealing with broken symmetries in the DFT of self-bound nuclei is a topic of active 
investigation, the best way forward is not clear. 



5.3 Single-particle energies 

The only quantities obtained from Kohn-Sham DFT for ground states that are guaranteed to be the 
same as those obtained from many-body wave function calculations (at an equivalent level of approxi- 
mation) are the total energies and the densities associated with the functionalso Of course, measurable 
quantities that can be expressed as the differences of ground state energies are also reliable. For nuclear 
applications, one would like to establish a robust connection between KS eigenvalues and nuclear single- 
particle energies, but it is often observed that Kohn-Sham eigenvalues, except at the Fermi surface, are 
not physical. On the other hand, Bartlett et al. claim that with a sufficiently rich Coulomb DFT func- 
tional, the KS orbital eigenvalues can be good approximations to removal (ionization) energies [58| [59]. 

We can understand how improved approximations can happen by considering how the full self- 
consistent one-particle Green's function G, whose spectral density is physical, can be expressed in terms 
of the KS Green's function Gks- We saw this earlier in the form of the Sham-Schliiter equation, Eq. fl65p : 
here we consider an alternative functional integral derivation and a diagrammatic representation to make 
some additional points. We add a non-local source ^{x',x) coupled to il){x)'il)'^{x') [we're in Minkowski 
space now with x = (x, t) so that we get real-time Green's functions]: 



t pijd'^x [C + J{x)ip''{x)iP{x)+Jd-^x'^{x)^{x,x')tP^x')] 



(220) 



Writing r[p, ^] = ro[p, ^] + rint[p, ^] and using functional properties of Legendre transforms, 
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which is represented diagrammatically in Fig. [13] |134j . (Note that these are the reducible self-energies 
here; so this is actually a rearranged Dyson-like equation that is equivalent to Eq. fp^ .) 
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A Gks 
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Figure 13: Full Green's function G in terms of the Kohn-Sham Green's function Gks- 

The Green's functions G and Gks yield the same density by construction; that is, the Kohn-Sham 
density pks(x) = —ii'G'^^{x,x~^) equals the full density p(x) = —ih'G{x,x~^). Here is a simple diagram- 
matic demonstration (the double line is minus the inverse of a single particle-hole ring) : 




^''Note that these densities are not necessarily obscrvables; e.g., the nuclear proton density is related to the measured 
charge density by a model dependent prescription, although this model dependence is generally considered to be negligible. 
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Figure 14: Schematic momentum occupation number n{k) for mean-field (Hartree-Fock) and with 
correlations. 

But other single-particle properties (e.g., the spectrum) are generally different, because the last two 
diagrams in Fig.[T3]will not cancel exactly. However, we can make them cancel more closely by redefining 
the Kohn-Sham system as described in Section [321 In the effective action approach, we find that adding 
sources does exactly this. For example, in the dilute Fermi gas EFT from Section the single-particle 
spectrum from a functional with p and r densities was shown to be closer to the exact spectrum (e.g., 
at the HF level) than a functional based on p alone (see Ref. [112j for details). This is consistent with 
the Bartlett et al. claim. More generally, we can use the Kohn-Sham basis in constructing the full G. 

One final note regarding Green's functions. The comparison of Kohn-Sham DFT and "mean-field" 
models often leads to misunderstandings, as when considering "occupation numbers", because of a 
confusion between G and Gks- Figure [H] suggests that occupation numbers are equal to or 1 if and 
only if correlations are not included. The Kohn-Sham propagator always has a "mean-field" structure, 
which means that (in the absence of pairing) the Kohn-Sham occupation numbers in the normal state 
are always or 1. But correlations are certainly included in T[p]\ (In principle, all correlations can 
be included; in practice, certain types like long-range particle-hole correlations may be largely omitted 
because of limitations of the functional.) Further, n(k) = (aj^Oj^) is resolution dependent (not an 
observable!); the operator related to experiment is more complicated. Additional discussion on these 
issues can be found in |172] . 

5.4 Improving empirical EDF's 

The technology for calculating with phenomenological energy density functionals such as those of the 
Skyrme form is already well developed and still improving. For example, machinery exists to calculate 
the entire mass table using a Skyrme form of the energy functional in a single day (examples of current 
capabilities are given in Refs. |106l [HI 1173] ). These calculations are quite successful, achieving root- 
mean-squared errors better than 2 MeV for the measured nuclides [H]. (See Ref. [35] for a state-of- 
the-art Gogny EDF.) Figures [TSl and [iGl show examples of both the successes and limitations of current 
Skyrme functionals. In particular, the trend is toward good reproduction of experimental data where it 
exists, but extrapolations toward the driplines where there is no data become uncertain and sensitive 
to poorly constrained parts of the functional. Improving the reliability of extrapolation is, of course, a 
prime motivation for ab initio DFT. 

There are substantial and diverse efforts in progress to improve the functionals in use. Most of these, 
even if not directly motivated by ties to microscopic ab initio input, will make such connections more 
likely. In particular these efforts include: 

• Generalizing the Skyrme functional. By making the functional a more complete representation of 
the possible physics, one gets closer to the systematics and model independence characteristic of 
effective field theory. Work in this direction includes adding tensor interactions |175l I113[ I176j . 
time-odd components |177[ 1178] . and gradient and density corrections [178j . 
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Figure 15: Deviations from experiment of Skyrme-Hartree-Fock EDF predictions for ground-state en- 
ergies by two functionals |174j . 
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Figure 16: Extrapolation of predictions from Skyrme-Hartree-Fock functionals and mass formulas for 
two-neutron separation energies in tin isotopes 



• Correlation studies. By studying systematically the errors from empirical functionals, clues for 
the important microscopic physics are unearthed. One example is the study of odd-even mass 
differences |173l 1179] and other examples are described in Refs. [HZl ESI I18U] . 

• Non-empirical approaches to pairing. At present these are hybrid calculations using a Skyrme EDF 
for the particle-hole part but a separable low-momentum potential at leading order to describe 
pairing [H] . The predictions of pairing gaps are remarkably consistent with those extracted from 
experiments (from three-point formulas using energies of adjacent nuclei as input) [Tl[|15[|16]. A 
systematic investigation of theoretical corrections is underway. 

• New constraints from ab initio theory. As noted earlier, microscopic many-body calculations of 
symmetric nuclear matter are less valuable than constraints from stable nuclei. However, the 
isovector parts of the empirical functionals are much less constrained by data and this is where 
the microscopic calculations should be most reliable. For example, Monte Carlo calculations of 
energies and densities for neutrons in traps (generally with AFMC). The trap serves as a variable 
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Vextlx), SO these calculations become benchmarks for testing or fitting empirical energy functionals. 
Note that constraints can be significantly more than realized by comparison to uniform neutron 
matter. Such a program is being carried out with GFMC/AFMC and NCSM calculations as part 
of the UNEDF project [9]. 

• Development of EDF without explicit underlying Hamiltonians. The strategy and distinction of 
EDF from DFT are described in Ref. [12], as well as challenges [13]. The insights gained and 
techniques being developed [121 SHI SS] will carry over to ab initio DFT. 

• Long-range correlations. With current EDF's, the incorporation of long-range correlation effects 
is handled separately from the functional calculations [HI SI]. The outlook for directly including 
this physics in an ab initio approach with OEP is unclear, but there are precedents in Coulomb 
OEP (e.g., the Coulomb ring sum describing the high-density expansion of the electron gas can 
be included). 

• Merging a variant of the density matrix expansion (DME) and empirical functionals. The com- 
parison of the DME-I and DME-II curves in Fig. [7] gives us an estimate of the truncation error 
in the expansion applied to the NNN terms because these prescriptions differ in the contributions 
of higher-order terms in the expansion. Indeed, it has been verified that suppressing these terms 
by hand brings the predictions for the A and B coefficients into agreement ^96]. The qualitative 
difference for the NNN-only contribution to B is large, but the actual coefficient itself is small, 
so this should not be alarming. However, because the combination of A and B and the kinetic 
energy to obtain the nuclear matter energy per particle involves strong cancellations, the spread 
in Fig. [7] is large on the scale of nuclear binding energies. 

These differences motivate a generalization of the Negele-Vautherin DME following the discussion 
in Ref. [181] . In this approach, the expansion of the scalar density matrix takes the factorized 
form 

p(R+ |,R- |) = ^n„(M(a(R)) , (222) 

n 

where 

(a(R)) = {p(R), r(R), VV(R), ■ ■ ■ } , (223) 

and kp is a momentum scale typically taken to be /cf(R) as in Eq. (11201) . Similar expansions are 
made for the other components of the density matrix. Input from finite nuclei can be used to 
determine the n„ functions, which can be viewed as general resummations of the DME expansion. 
A program to merge these developments with empirical Skyrme functionals is underway ^115j . 

Almost all of the major developments are only recently begun, so many new results can be expected in 
the coming years. 

6 Summary and outlook 

The quest for an understanding of all nuclei based on microscopic few-body interactions among nucleons 
is being attacked from multiple directions. In this review, we have considered ab initio density functional 
theory for nuclei, defined as DFT based on a realistic nuclear Hamiltonian (meaning one that accurately 
reproduces few-body data) and orbitals that satisfy equations with local Kohn-Sham potentials. This is 
an ambitious undertaking and one which will not be realized without significant further developments in 
nuclear structure calculations. Nevertheless, the prospects are good for making microscopic connections 
to empirical functionals in the short term and building on them steadily. 
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We base this optimism in part on the successful advances in ab initio DFT for Coulomb-based 
systems by quantum chemists and condensed matter physicists. We believe nuclear physicists will 
profit from greater attention paid to Coulomb DFT developments. At the same time, we recognize 
that these advances have taken many years to realize and may not be readily transferred to the nuclear 
domain. Indeed, we have argued that strong analogies between Coulomb and nuclear systems are really 
only apparent when one considers low-momentum interactions, so that the correlations induced by more 
conventional (local or almost local) interactions are tamed from the outset. As one strives for greater 
accuracy, the differences between atoms or molecules and nuclei are likely to become more significant. 

Perdew and collaborators describe for Coulomb DFT a "Jacob's Ladder" of increasingly sophisticated 
"nonempirical" functionals stretching toward the heaven of chemical accuracy. In Ref. [12], five rungs 
of the ladder are identified (see also |182] ): 

1. The local density approximation (LDA) with both spin densities P|(x) and pj,(x) as ingredients 
(usually called LSDA). 

2. The generahzed gradient approximation (GGA), which adds dependence on Vp|(x) and Vp|(x) 
to the LSDA functional. 

3. The meta-GGA adds dependence on (some subset of) V^p|(x), V^p|(x), T|(x), and Tj(x). Note 
that the kinetic energy densities are actually non-local functionals of the ordinary density, although 
only semi-local functionals of the oocupied KS orbitals. 

4. Hyper-GGA, which includes the exact exchange energy density calculated with the (occupied) 
orbitals. 

5. Full orbital-based DFT, which in addition to exact exchange uses unoccupied orbitals. For ex- 
ample, this might include the random phase approximation (RPA) with Kohn-Sham orbitals to 
address long-range correlations. 

An important aspect of the overall approach is "constraint satisfaction" as opposed to data fitting, which 
is why the functionals are called nonempirical [182j . Climbing this ladder has been a decades- long effort 
by quantum chemists with much development in progress on the last rung. 

We can imagine a corresponding ladder for nuclear DFT climbing toward a universal energy density 
functional tied to microscopic nuclear Hamiltonians (and ultimately to QCD) that predicts known 
nuclear properties more accurately than currently possible and robustly predicts unknown properties 
with plausible theoretical error estimates. The rungs might look like this: 

1. Present-day Skyrme EDF's, which are mostly empirical (i.e., fit to properties of selected medium- 
to- heavy nuclei), such as described in Ref. [TTj. 

2. Generalized Skyrme interactions, with additional gradient and density dependences, with new 
constraints from microscopic theory (e.g., neutron drops). 

3. Functionals that merge the long-range parts of the microscopic NN and NNN interactions pre- 
scribed by chiral EFT, converted to semi-local form with a variant of the density matrix expansion 
(DME), with a Skyrme functional. The short-range parameters should be refit to nuclear proper- 
ties and theoretical constraints. 

4. A complete functional based on a variant of DME applied to a low-momentum potential that is 
evolved from chiral EFT NN and NNN interactions. 

5. Full orbital-based DFT based on low- momentum interactions. 
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This ladder is tied to our restriction to local Kohn-Sham potentials and so naturally builds on Skyrme 
EDF's. An alternative ladder could build on non-local potentials (e.g., on Gogny-type EDF's) that arise 
from derivatives with respect to density matrices rather than densities; see Refs. [T^ [TB] for discussions 
along these lines. Just as in the Coulomb ladder, one hopes for monotonic improvements at each rung, 
but this may not always be the case. (Of course, rungs will shift or more will be added as progress is 
made.) 

Rather than climb slowly over the course of decades, nuclear physicists are trying to catch up rapidly 
by attacking all of the higher rungs in parallel [Hlinilin]- As described in Section [531 there are extensive 
ongoing efforts to explore the second rung, including more general density dependence and higher-order 
gradients [9|,[T0], while projects on both rungs three and four have recently been started as part of the 
UNEDF project |8]. The last rung is in the exploratory stage with tests started in one dimension with 
nuclear-like interactions for both self-bound and trapped systems to explore the accuracy of the KLI 
(and related) approximations compared to a full OEP treatment of such systems and to establish the 
limitations of semi-local expansions such as the DME. 

While the last two rungs can be purely predictive if only few-body data is used to determine the 
input Hamiltonian, it is likely that fine-tuning to heavier nuclei will be needed to reach the accuracy 
goals. This is because even the best ab initio nuclear structure calculations at present only achieve 
about 1% accuracy in ground-state energies for a small fraction of the table of nuclides, while nuclear 
EDF's such as Skyrme functionals easily surpass this over most of the table. Indeed, the comparatively 
high accuracy of nuclear EDF may imply that a different organization of the problem (e.g., based on a 
finite-density effective field theory of nuclei) may be needed. 

There are many important open questions to be addressed in the course of these projects (see also 
Ref. [IS])- Among them are those highlighted in Section [S]on symmetry breaking and restoration, single- 
particle levels, and pairing; these will be relevant for all rungs of the ladder. At the top of the ladder 
we can ask: Which of the problems encountered in Coulomb DFT that motivate orbital-based func- 
tionals have analogs in nuclear DFT? E.g., what is the impact of self-interaction and missing derivative 
discontinuities? Another class of open questions concerns the input Hamiltonian: For low-momentum 
potentials, how accurate will many-body perturbation theory be? What corrections/summations will 
be needed to reach the desired accuracy for nuclear DFT? Are four-body forces necessary? The list of 
problems yet to solve might be intimidating, especially given the many decades already spent attacking 
nuclear structure, but the coherence of the current effort and the dramatic advances of just the last few 
years give hope that they can be overcome. 

We have taken a broad view at what ab initio density functional theory could be but our discussion 
has been far from an exhaustive treatment. By focusing our discussion on orbital-based DFT for 
nonrelativistic Hamiltonians we have excluded various alternative paths to ab initio DFT. Here are 
some of the other approaches that are relevant to the wider nuclear DFT effort: 

• The perturbative in-medium results from low-momentum potentials suggest that an alternative 
EFT power counting may be appropriate at nuclear matter densities. Kaiser and collaborators 
have proposed a perturbative chiral EFT approach to nuclear matter and then to finite nuclei 
through an energy functional |183[ I109| IllH [TMt 1185] (see also [186] ). They consider Lagrangians 
both for nucleons and pions and for nucleons, pions, and A's, and fit parameters to nuclear 
saturation properties. They construct a loop expansion for the nuclear matter energy per particle, 
which leads to an energy expansion of the form 

oo 

Eikp) = ^2^^^ fnikp/m^, A/m^) , [A = - Mn ^ 300 MeV] (224) 

n=2 

where each /„ is determined from a finite number of in-medium Feynman diagrams, which incorpo- 
rate the long-distance physics. All powers of k-p/rriT, and A/m-,, are kept in the /„'s because these 
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ratios are not small quantities [187j . A semi- quantitative description of nuclear matter is found 
even with just the lowest two terms without A's and adding A's brings uniform improvement 
(e.g., in the neutron matter equation of state). There are open questions about power counting 
and convergence, but many promising avenues to pursue. By applying the DME in momentum 
space to this expansion they derive a Skyrme-like energy functional for nuclei, which has also been 
merged with a treatment of strong scalar and vector mean fields. |109[ IllOl lllll 11841 1185j . 

The Superfluid Local Density Approximation (SLDA) and extensions developed by Bulgac and 
collaborators [1491 1188^ 11891 1190] have been applied with great success to systems of trapped cold 
fermionic atoms with large scattering lengths. The SLDA is an extension of DFT to superfluid 
system that uses a semi-local energy density with parameters determined by matching to numeri- 
cally accurate Monte Carlo simulations of uniform systems. Extensive further development of the 
SLDA for nuclear systems is in progress. 

We have noted that the Kohn-Sham DFT emphasis on locality in coordinate space is not so natural 
in many-body formulations. In particular, the nuclear shell model is naturally considered in the 
space of orbitals. Papenbrock has applied the Legendre transform formulation of DFT to energy 
functionals based on shell-model occupation numbers. That is, a functional of the density now 
becomes a function of the occupation numbers in the model space. He and his collaborators have 
shown the usefulness of this in solvable models (the pairing Hamiltonian }191] and the three-level 
Lipkin model [56]). The functional is orbital-based and thus non-local in density. For applications 
to nuclei, the idea is to generalize the Dufio-Zuker mass formula |1921 1193] to a functional of 
occupation numbers, with values determined from the minimization of the functional. With 10 
parameters, an RMS value of deviation for masses of about 2000 nuclei is just over IMeV. 

Schwenk and Polonyi have proposed an alternative approach to ab initio but orbit-free DFT 
(i.e., not Kohn-Sham) using a clever renormalization group (RG) evolution [132] . The idea is 
to introduce an effective action for a nucleus with a low-momentum interaction included with 
a multiplicative factor A and a confining background potential (e.g., a harmonic oscillator trap) 
with a factor (1 — A). As A flows from to 1, the background potential is turned off and the 
interactions, with associated many-body correlations, are turned on. This evolution is dictated by 
an RG equation in A. Work is in progress to implement this approach for a realistic nucleus |194j . 

There is a large body of work on covariant approaches to nuclear DFT, which cannot be adequately 
summarized here. Fortunately, there are many reviews that highlight DFT connections, including 
Refs. |195[ 11961 1197[ 11471 1198[ I199j . to supplement our brief discussion. Covariant DFT invokes a 
different organization of the nuclear many-body problem that is particularly compelling because 
of the coupling of spin-orbit and central components dictated by relativity. In particular, the 
characteristic feature of relativistic approaches to nuclei is large isoscalar Lorentz scalar and 
vector mean fields, several hundred MeV in magnitude at nuclear densities, which cancel to 
provide nuclear binding but add to produce spin-orbit splittings. 

Relativistic mean-field models were originally motivated as the Hartree approximation to a more 
complete theory, but a more recent interpretation of the largely empirical functionals is as rela- 
tivistic versions of Kohn-Sham DFT [196] . The mean fields serve as local Kohn-Sham potentials 
in Dirac equations for the orbitals. This picture is consistent with a microscopic derivation within 
an effective action formulation as in Section HI with sources coupled to each of the relativistic 
densities or currents or to corresponding meson-like auxiliary fields [147] (this is in contrast to 
implementations of relativistic electronic DFT for heavy atoms, which include only the vector 
four-current |3]). Application of a loop expansion with proper renormalization gives the possi- 
bility of a systematic microscopic expansion [2001 1201] . Other ongoing efforts to provide more 
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ab initio aspects including a connection to free-space NN scattering derive relativistic EDF with 
density dependent meson-baryon vertex functionals, derived from the Dirac-Brueckner theory (see 
Refs. 12021 EOSl 1203 and references cited therein). 

These more microscopic approaches do not at present have the close connection to few-body data 
and to ab initio structure calculations of light nuclei that is possible with the non-relativistic 
approaches we have considered. However, the wide phenomenological successes and complemen- 
tary nature of covariant nuclear DFT motivate further work toward incorporating additional 
microscopic constraints. In the process, most of the challenges confronting nonrelativistic DFT 
discussed above are being attacked in parallel for covariant DFT. This includes the issues of real- 
istic pairing interactions |205i 12061 1207] . coupling to particle-hole vibrations [208j . and symmetry 
breaking/restoration [203 12101 1211] • 

Each of these avenues has the potential to make a significant impact on nuclear DFT. 
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